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Summary. The difference methods for the numerical integration of ordinary differ- 
ential equations of the first order are discussed by using operator calculus and symbolic 
expansions. A new straightforward central difference method is developed, which is 
based on a formula closely associated with Simpson’s rule. The main features of the 
method are that, for each step of integration, the largest unknown term is determined 
by an algebraic equation and that the remaining difference correction is extremely 
small. The method can directly be applied even to systems of the first order with one- 
point houndary conditions. A numerical example is given. 

1. Introduction. Different methods have been given in the literature for the numerical 
integration of the ordinary differential equation of the first order 


= F(a, ») (1) 


with the initial condition y = yo for x = 2 . The purpose of this paper is to discuss 
only the methods based on the theory of finite differences, since they seem to combine 
simplicity with a large general applicability. The fundamental formulas are here of 
two types: central difference (CD) formulas and backward difference (BD) formulas. 
The central difference formulas, due to Gauss, are characterized by an extremely rapid 
convergence or semi-convergence and a small error-term, but in using them each step 
of the integration demands estimation and iteration. The use of estimations is avoided 
altogether in the integrations by means of the backward difference formulas, of which 
the first is due to Adams-Bashforth, but instead are these formulas slowly convergent 
and have relatively large error-terms. Both these methods have difficulties in knowing 
how to start, and in general it has been recommended to start the integration inde- 
pendently by means of a Taylor-series expansion. 

By using operator calculus and symbolic expansions, the connection between the 
CD-formulas and the BD-formulas will here be investigated in greater detail. Utilizing 
the experience found in this way, we will show that it is possible to construct integration 
methods which combine the straightforwardness of the BD-methods with the simplicity 
and rapid convergence of the CD-methods. Even the starting problem will be simply 


solved. 


*Received July 25, 1951. 
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2. Difference operators. The extrapolation principle. Let h be the interval, let E be 
the step-operator defined by Ef(x) = f(x + h), and let f, mean f(a + nh). We will 
then introduce the operators A, V, 6, and uw for the formation of forward, backward, 
and central differences, and mean values, respectively, by 


A=E-1, V=1-E", (2) 
é=f’-E™, p= (E'" + E”)/2. (3) 

From (3) one obtains directly 
w=1+8/4, wl+ 0/4"? =1. (4) 


which relations often can be used in transforming CD-formulas into a suitable form. 
In the following we will use operator calculus and symbolic expansions’ (Sheppard 1899; 
see also Michel 1946, Bickley 1948). If D = d/dz is the differentiation operator, Taylor’s 
series gives EH = exp (AD), and from (3) we then get 


hD = 2sinh™ 6/2, (5) 


which is the basic CD-formula for numerical derivation and integration.’ 

Among computers it is now a well-known fact that, in using pure CD-formulas or 
BD-formulas for some purpose, it is often impossible to utilize all the function values 
in a given material, which is illustrated by the first two figures below: 


NSAAASAAASSASS AAS 











Wasa 


BD-formula CD-formula CD-formula with some CD extra- 
polated by means of available BD 

(or mixed CD-BD-formula) 
The triangle indicates a function given numerically in equidistant points and its difference scheme. The 
full line shows the differences involved in a difference formula of a certain type; the part of the function 
values taken into account in this way is shaded. The dots in the last figure indicate CD extrapolated by 


means of available BD. 


Bickley and Miller (1942) pointed out that there exist an infinite series of “mixed” 
CD-BD-formulas by means of which a given material could be taken into full account, 
and they have worked out extensive tables for the numerical derivation. As far as we 


1The formulas obtained in this way have only symbolic character, but they can be rigorously derived 
and their remainder-term can be determined by starting from Newton’s interpolation formula for un- 
equal interval with Cauchy’s remainder (see Nielsen 1908). 
2For tables of coefficients, see e.g. Salzer (1943, 1944, and 1945). 
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know, corresponding tables for the numerical integration have not yet been published, 
nor is it necessary. 
Here we will proceed in another way. From (2) we obtain 


E=14+VE=14+ 7404-4 V+ VE 


(6) 
2 3 -1 
=14+V4+V4+V7V+°-:-=(-YV) , 
and according to (2) and (3) we then get the extrapolation formuias 
oF, = "1 — 9) "Fe, = Dall 7 F,,, , 
k=0 
5 mel _= * eae | ies Vv)’ m » - a 
k=0 
(7) 
el , 2m/<« r-m— —m Qm+ek 
ud" Pease = V2 — VL — FV) Face = Dy OT Fase 5 
k=0 
‘ 2m ¢ r—-m— 7 —m 2m +k 
a"? = FP" = TN = OF he * 2 eae 
k=0 








.| 
‘Oo 1 23 4 \jo 1 2 3 4 
0 1 O 0 0 0 ei2 1 1 l i 
~1/11 11 1 -1123 4 5 6 
2/12 3 4 5 2/25 9 14 2 
a i: 

3/13 6 10 15 - -312 7 16 30 50 
4/1 4 10 20 35 - -~412 9 25 55 105 





satisfying elementary recurrence relations. The formulas (6)-(7) form together an 
“extrapolation principle”, by means of which CD-quantities in the horizontal lines n 
and (n + 3) can be expressed in terms of the BD in the backward line (n + r). The 
series in (7) are only formally infinite; in the practical use they are always interrupted 
after a difference of the finite order p, corresponding to an extrapolation of the function 
under consideration by a polynomial of the degree p. 

Instead of using the rather complicated “mixed” CD-BD formula of the Bickley- 
Miller type, we can now take a given material into full account simply by using a pure 
CD-formula and by extrapolating as many additional CD as possible by means of the given 
BD in the last backward line available. Since the coefficients in (7) are all integers, the 
extrapolations can be rapidly carried out on ordinary desk machines. The result will be 








100 PER-OLOV LOWDIN [Vol. X, No. 2 


the same as obtained by using a ‘‘mixed”’ formula, and we have to remember that the 
error in our case is given by the sum of the remainder in the pure CD-formula and the 


remainders in (7 
3. Gauss’s formulas and the classical backward difference methods. According to 
(4) and (5) the first orden equation y’ = F(z, y) is equivalent with the CD-equation 


hF = p(] > & 4) a sinh ; 6 2-y 


(S) 
; res Le 
_ e a. wi oS OS 07 SP 62 
ms gh eet a94 °F 1494 °F ?, 
from whit h y can be solved by the inverse relation 
h7'y = p(1 + 6°/4 ‘2 sinh’ 6/2!~'-F 
= (9) 
pai ne 37 191 ssp 4 2497 a 
uu Oo ees” L or a JS Jia Ey a : oc ™ 
yet  aae 60480 * ‘ 3628800 % 
Here the symbol 6 'F means the first sum of F, defined by 6(6 'F) = F; the “summation 


constants” included in this and other sums will here be determined by the condition 
that the integration formulas under consideration should be valid even in the starting 
point a z, . Formula (9) forms the basis of a method of numerical integration” due 
to Gauss, which was first published by Encke (1837) and which has frequently been 
used in the celestial mechanics (see e.g. Charlier 1907). 

Here we will study a slightly modified form. By applying the operator 6/’~ on both 
members of (9) we get 

{ 

he by sal onl : 3! Basie F =o OP sie we eFectt +e* EO) 
a form used for instance by Hartree (1928). This is nothing but the (rapezordal rule 
with difference correction, and we note that, in integrating a given fix integrand F = 
F(z), tl 
princi 
show that it is possible to avoid the use of estimations, characteristic for the earlier 


1e whole material can be taken into full account by using the “extrapolation 


le’ developed in §2. In integrating the differential equation (1). we will now 


CD-methods, by using the same principle. 
Let us assume that we have started the integration in some way and that we have 


computed y, F = F(x, y), and the difference scheme for F accurately up to the point n 
given by x = a + wh. Since we know the BD of F in the backward line n, we can now 
extrapolate the CD in the horizontal line (n + 4) in formula (10) according to (7), 
and in this way we find approximate values of y,,, and F,,,, . Repeating the process 
we obtain approximate values of y,42 , Piso. Yass » Pass , +++ - Going back to the point 
n + 1), we can then improve the accuracy of the solution by the repeated use of 


formula (10) and the extrapolation formulas (7), utilizing the CD found on later stages 
of the first approximate calculation and extrapolating some additional CD by means 
of the last BD-line available. By iterations it is in this way possible to obtain a pure 


CD-result in a straightforward manner. 


’The remainder was first derived by Nielsen (1908); see also Steffensen (1924) and Nystrém (1926). 
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The connection between this method and the classical BD-methods is perhaps of 
some interest. Putting the expressions for the CD in (10), extrapolated according to (7) 
by means of the BD-lines n, n + 1, n + 2, --- , into formula (10) we get 


=} 7 I ’ 5 on 3 eal 
ho'y — hy, = Prat 5 VPn t+ ig VP t+ gV 
+h vv, += vr, + 
720 . 
| | as) 
I 1-9 Vi i 19 V Pasi — og VF , 
i i 
709 V Paes 160 ¥ fat - ' 


which is nothing but the BD-formulas given first in another way by Adams-Bashforth 


1883). The accuracy of the first formula is low,* corresponding to our first approximation 
above, and it is usually recommended to improve the accuracy by using one or more 
additional BD-formula, which may formally be obtained from the ground formula by 
letting the operator 1 = (1 — Y)E work repeatedly on its right-hand member.’ A 


comparison between the formulas (11) and the simple formula (10) shows immediately 
that, from the point of view of the computer, it is considerably simpler to use the single 
CD-formula (10) together with the extrapolation formulas (7) than the corresponding 
classical BD-methods. 

Let us now also consider formulas for taking a double step® from (m — 1) to (n + 1) 
by means of the values associated with n. Letting the operator ué work on both members 


of (9) and using (4), we get in the point n 


isis ys lap l wp 4 1 yop 23 asp 
l iY, = hh, - Pa = 1 —— fhe oe ry 
M rg6°*=~ ia9°"** 1519 °"* saeg00"** 


+ -::, (12) 
which is nothing but Simpson’s rule with difference correction.’ By combining (12) with 
the extrapolation formulas (7) and using iterations, we can again construct a straight- 
forward CD-method in the same way as described in connection with formula (10) for 
taking a single step. 

The connection between this method and the classical BD-methods is easily seen. 
By putting the expressions for the CD in (12), extrapolated according to (7) by means 


Che remainder was first derived by Nielsen (1908); the propagation of errors has been investigated 
by v. Mises (1930) and by Schulz (1932). The method has been further developed in Russian memoirs 
by Krvloff nodification has been given by Falkner (1936). 

See also the note bv Stohler (1948). 

’The first method of this type was developed by Richardson (1911, 1927), who used approximate 
CD-expansions and afterwards improved the solutions by a certain process called the ‘‘h?-extrapolation” ; 
mpare also Dunean (1948). 


'The simple Simpson’s rule has earlier been used for numerical integration by Bickley (1932) and by 
Collatz and Zurmiihl (1942a). 
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of the BD-lines n, n + 1, n + 2, «++ , into formula (12) we obtain 


ho"Ynar — Ro"'Yn-1 


| - 29 14 ‘ 
OF. +=VF, +577. + — VP Fr. 4 - 
3 b 3 VE. + 99 Vn + a5 VE + (13) 
ae ee | ee) 
~aens vo a ale iti 90 did = 


The first formula (13) was derived by Nystrém (1926), and later a slightly modified 
form was given by Lindeléf (1939); his correction term can here be obtained by taking 
the remainder terms in (7) into account.” The second formula (13) has been treated by 
Levy-Baggot (1934) and by Sibagaki (1936). Again we find that it is simpler to use the 
single CD-formula (12) together with the extrapolation formulas (7) than the corre- 
sponding classical BD-methods. 

The main conclusion in this section is therefore that, from the computer’s point of 
view, it is considerably simpler to use a fundamental CD-formula in combination with 
the extrapolation principle (7) than any of the classical BD-methods. However, if high 
accuracy is desired, the numerical integrations based on (10) or on (12) are rather 
laborious to carry out, since each step of the integration demands iterations of the whole 
difference scheme. In the next section we will therefore try to modify the basic CD- 
formula in order to find a still simpler and more straightforward integration method. 

4. A new central difference method. 1) 7HE BASIC FORM ULA. Let us first try 
to derive a new formula of the Gaussian type (9), but having a much smaller difference 


9 


correction. According to (3) and (12) we have 


; ] ] 23 
at a re a 5h o— UF Pe Al me 67 __ a oe ae = (14) 
hp by = F 4 3 (H I VF) 180 5F + 1512 5} 06800 ° | + 14) 


Letting the operator (ué)~* work on both members of this equation, using (4) in treating 
the difference correction, and transforming the F-term to the left member, we get the 
integration formula 


hy — 3 F(a, y) = (ud) "{F — ; VF} 


(15) 
31 : 557 


] —— : an _ 5074 
oF + F — 907200 


180 © 15120" ° ee 


which forms the basis for our method. The sum in the right member is dominating and 
will be called the main term: 


se Wiad : 
M = (u5) F - = VF?. (16) 
38) 

The difference correction is extremely small and will be denoted by y. The quantities 


involved in the computations are suitably arranged in two or three tables: a main table 


8See formula (6), first line. 
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recording the quantities z, y, (M + y), y, M, and F, a difference table giving the differ- 
ences of F’, and—for non-linear equations (1)—an iteration table. As to the difference 
correction the integration will be divided in three. parts: the start, the marching process, 
and the aftercorrection. 

Il) THE MARCHING PROCESS. At first we will assume that the integration has 
been performed up to the point n with 2 = 2) + nh. In order to carry out the next step 


of the integration, we will then write eq. (15) under the form 


Mas M,_, + 2F, — VF, (17’) 
a + - oI oF. —- a ar de oes (17’’) 
| atte 130 4 ° “=+ T 75199 4 9 *=*1 ~ 997900 4 ° Fs , - 
1 . 
lh Y , 9 F, , = (M f VY) n+1 i (17’”’) 


First the new main term ,,, is computed from the quantities in the previous part of 
the main table according to eq. (17), which is obtained by letting 2ué work on both 
members of (16). Then we will extrapolate’ the CD occurring in the difference correction 
Yn+1 by means of the known BD in the backward line n by using formulas (7): 


Qu 6°F,., = 2V°F, + 5V'F, + 9V°F, + 14V°F, + --:, 
(18) 
Qu 6F,., = 2V°F, + 77°F, + 16V'F, + 30V°F, + ---, 


Finally we will determine the values of y,,, and F,,,, by solving the algebraic equat ion’? 
(17’’). For linear equations (1) the solution is obvious, and for non-linear equations 
we will recommend the use of an iteration process, which will be discussed later (VI). 
The quantities y,,, and F,,,, are written down on their places in the main table, a new 
backward line (n + 1) is calculated in the difference scheme for F, and then the process 
can be repeated for the point (n + 2) with « = 2) + (n + 2)h, ete. This integration 
process is quite straightforward, and it gives a preliminary solution y with a high 
accuracy, due to the smallness of the difference correction. 

il) THE AFTERCORRECTION. The preliminary solution obtained in the march- 
ing process has an accuracy which is determined by the “mixed’’ CD-BD character of 
the integration formulas (17) and (18). However, when the marching process has been 
finished for the whole range of integration under consideration (including some addi- 
tional points), it is possible to improve the accuracy of y by utilizing the actual values 


°A necessary condition for the usefulness of such extrapolation is, of course, that the higher differences 
| of F converge rapidly against zero; this can be obtained by choosing the interval h sufficiently small. In 
general we will calculate the differences of F only up to the order where the irregularities, due to the effect 
of rounding-off errors, become to be of about the same magnitude as the differences themselves. 

Tn agreement with Hartree (1949, p. 235) the term “algebraic” is here meant only as an antithesis 
to the term “differential.” 
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of the CD found on later stages of the calculation. For this purpose we will introduce 


the correction 


corr, = y, , actual — y,,, extrapolated 
= (19) 
3 ry o om 05 r7 
= ——— fy OF =~ OF unl + ZH Sg a SY OY 
io . TY 15120 ' sti " 


We note that, in the main table, the quantities y and WV are usually recorded with the 
same number of significant figures, whereas it is often sufficient to record F with at 
least one figure less'* than M. The correction (19) will, in general, influence the value 
of y, but the idea of our method is now that the interval h should be chosen so small, 
that the change in y should not influence the recorded figures of F. In this case the 
columns for F and .V/ in the main table will be unchanged by the aftercorrection, and, 
according to (17’’’), the final solution can be computed from the simple formula 

Yn,final — Ynys nary + Acorr, . (20) 
In this way it is possible tO carry out the transition from a mixed CD-BD-result into 
a pure CD-result by a single iteration of a very simple type (compare §3). 


During the marching process it is suitable to check the extrapolations (18), e.g. in 


every fifth point, by comparing the extrapolated value of 2ué°F, in the point n with 
the actual value, found when the integration has reached the point (n + 2). If the value 
of h seems to begin to become too large, it is then necessary to change the length of the 
interval by some of the well-known processes of subtabulation. 


IV) THE START. Even the start will here be treated by means of our basic formula 


(15) by using an auxiliary formula’ and a method of successive approximations. Letting 
the operator A = ywé + 6/2 work on both members of (15) and using (4), we get 
a a ee or a 
AM =F+ i 6F¥ = F+ a6 T av 6F, (21) 
8) © \ | 
which for x = 2, gives the expansions 
; = l a in ie ips 
B= iat, See, ~ oe 8, + S487, ~' 4 
6 24 YO 
(22) 
l ; 5) I 
I 77 6 hF3 + 97 oi 2] 0Yo T 


where the latter is derived from the former by using (5) and (17); the second form is 
the best for our purpose. The integration will now be started from the initial condition 


y = yo for x = 2, and by using the formulas 

| 

i, = hy =F, —4 
o 

(23) 

2 h 3 l 

Ms = Wy +2 + ER + oy tx hy t 

1("“ompare footnote 9 


~Compare Richardson and Gaunt (1927). 
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In a zero-order approximation, the quantity y is entirely neglected: y = 0. The quantities 
M, and W, found from (23) are written down on their places in the main table, and 
then the quantities y, , F,; , Ms, y., F2 , Ms, «++ ete. are computed successively by 
using (17’) and (17’’). In this way some approximate values of y and F, forwards and 
backwards'* from the starting point, are determined. 

Now the difference scheme for F provides a first-order approximation of y in the 
neighbourhood of the starting point. Using these values of y, we can then repeat the 
start until two consecutive approximations of y agree within the significant figures; the 
integration can then be continued as a marching process (II). Due to the extreme 
smallness of the difference correction, the first order approximation of y shows a sufficient 
accuracy in most problems. 

In problems where the functions involved in eq. (1) are specified only in one direction 
from the starting point, it may be necessary to extrapolate the CD in the range of 
x = 2 by means of forward differences, using a particularly small interval. 

The method of successive approximations described here can easily be extended to 
the whole range of integration,'* but, if the computations have to be carried out only 
by the aid of desk machines, it is certainly suitable to confine the method to the start. 

VY) THE AUTOMATIC ELIMINATION OF OSCILLATING ERRORS. A partic- 
ular phenomenon in the first approximation of the start will here be briefly mentioned. 
Due to the neglection of y» in (23), the approximate values of 7, and /, are affected 
with small errors of about the magnitude —y,) and +1.5y, . This implies that there is 
also an oscillating error in the first difference scheme for F, which is propagated to the 
higher differences with a steadily increasing magnitude, and the difference scheme can 
therefore have a rather irregular appearance. '” 

At first sight it seems impossible to determine y with a sufficient accuracy from a 
difference scheme which is disturbed by an oscillating error. However, it is easily proved 
that an error of the type (—1)"f, , where / is a polynomial in x of the degree p, can be 
entirely eliminated by the operator y’ 


pw *''(—1)"f,} = 0. (24) 
The operator » can be introduced in any difference formula by the repeated use of the 
unity operator 1 u (1 + 6 /4)~', and we note that the application of the refinement 
process 
, hid ir Di lean , 6 ae ox 
eon ape ——p OF +—_ae OF —- —Hw OF + —-+:. (25) 
| 16 64 


will every time diminish the degree of the oscillating error by two. 

In our case, i.e. the first approximation of the start, the magnitude of the error is 
very slowly varying, and we can therefore conclude that, due to the appearance of the 
mean-value operator uw in (17’’), the oscillating error will automatically be almost entirely 
eliminated in forming the difference correction y. In the second approximation the 
phenomenon has therefore disappeared. A numerical example may be found in Table IT. 

8When working backwards, we obtain the formulas needed by substituting —h instead of h in (15); 
we note that there are different values of M, for the two directions. 

“Compare the treatments of linear equations (1) by Hausmann and Schwarzschild (1947) and by 
Fox and Goodwin (1949). 

The existence of oscillating errors has been reported also by other authors (e.g. Fox and Goodwin 
1949) in connection with other methods. 
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VI) THE SOLUTION OF THE ALGEBRAIC EQUATION. The straightforward- 
ness of our integration method based on (15) is essentially depending on the fact that 
we have introduced a certain algebraic'® equation, by means of which the largest un- 
known central term F is determined for each step of integration. A similar idea has previ- 
ously been used by Noumerov in treating second order equations. The eq. (17) may 


be written in the form 
y V + cF (a, y), (26) 


where .V hA(M + y) and c = h/3. For linear equations (1) the solution of (26) is 
obvious, and we will therefore confine the following discussion to the non-linear case. 

A well-known analytical solution of (26) is given by Laplace’s formula, but, from 
the practical point of view, we will here instead recommend the use of an iteration 
process. Let us first consider the marching process (II), and let us assume that the 
integration has been performed up to the point n. According to (6) a zero-order approxi- 
mation F,,..; of F,,, ean then be found by extrapolation: 


Fe. =Fi+VFA.+ VFA+::- + VF, . (27) 


By putting this value of F,,,, into the right-hand member of (26), a first-order approxi- 
mation of 7,4; is obtained; then a new value of F,,,, is calculated ete., until two con- 
secutive approximations of y,,, agree within the significant figures. The solution is 
illustrated by 

/ V+cF’’; 9" N+cF"”; yo N+cF"”;.--- (28) 
where the upper index gives the order of approximation. This iteration process is of 
the first order,’’ and the condition for convergency is | cf, | < 1. Since the zero-order 
approximation of F is in general rather accurate,'* the solution is rapidly obtained. A 
numerical example may be found in Table IV. A check on the calculations in the differ- 
ence scheme. is provided by the relation 


oO Pcs Bice Fens (29) 


Let us then consider the start (IV). From the practical point of view we will here 


recommend the use of a second order process of a simple type. Since the process (28) is 


of the first order, the differences y” 2 y, °*: ete. form approximately 
a geometric series’ with the quotient ¢ = cf, . Summing this series, we get 
yr — yy" 
- 1 ‘ 
> y tly yl + io + 
t~*< 
(30) 
l } 
Y a ae 
Y ; - : 
y -2y" + y 


‘6Compare footnote 10. 

"Iteration processes of this type were first classified by Schréder (1870); see also Hartree (1949). A 
formula due to Theremin (1855) corresponds to a process of order 
This fact forms the basis of all the classical BD-methods. 
This is easily proved by putting y = y + 7»™, and expanding F'(z,y) into a power series in »™. 
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where the middle step for gq = cF§” corresponds to the Newton-Raphson formula.” 
We will instead use the formula given in the last step. The value y* of y found in this 


way can be checked by a new iteration y** = N + cF*. If a still higher accuracy is 
required, the whole process can be repeated by introducing the values of y*, y**, and 
yt** N + cF** in (30). This manner of proceeding corresponds to a well-known 


‘ . 1 
iteration process of the second order.” 


Putting (28) into (30), we get the final formula 
si TT ai iat 
7" “* saa (2)? 
Pp — 2" + F 





(31) 


which seems to be the best for our purpose. In using this formula, F“” is determined 


from the few terms available in (27), and then the final result can be obtained by means 
of only two computations of the function F = F(z, y);for a numerical example, see Table V. 

The results obtained in this section show that the rather laborious iterations of the 
whole difference scheme for each step of integration, characteristic for the classical CD- 
methods and for many of the BD-methods (§3), have here been reduced to the after- 
correction (20) and to the simple iteration processes described above for solving the 
algebraic equation (26). For linear equations (1) the simplifications are still more con- 
siderable. 

VII) CHECK OF THE SOLUTION. In all step-by-step methods based on the use 
of recurrence relations, it is necessary to have an accurate check on each stage of the 
calculation, since a mistake somewhere will vitiate the whole subsequent work. The 
most important check is here provided by the difference scheme for F, since a small 
error in the solution will give rise to a large irregularity in the higher differences of F; 
in this connection even eq. (29) may be of value. By letting the operator 2ué6 work on 


both members of (15), we get 


h 


Yn+i — Yau = 3 (Fast + 4F,, + F,-1) + 2hu On ’ (32) 


which relation may also be used for check purpose. A more independent check is provided 
by the differential equation itself in the form (8). 

VIII) NUMERICAL EXAMPLE. In order to illustrate how our integration 
method works in a practical case, we will here give some results from a treatment of 


the non-linear equation 


dy 2 ‘ 

oa eo of 33 

dx ie (33) 
with the initial condition y = 0.72 901 1133 --- for x = 0. The exact solution is here 


the logarithmic derivative of the Airy integral 
© log Ai(2) (34) 
y = — log/ 
y= 7, “San 


which has been carefully tabulated by Miller (1946). The start of the integration is 
illustrated in Table I, the elimination of oscillating errors in Table II, the marching 


2Compare Collatz and Zurmiih] (1942a). 
Aitken (1925), Samuelson (1945), and Hartree (1949). 
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process in Table III, and the solution of the algebraic equation in Tables IV and V. 
The calculation was carried out with nine figures with a speed of 12-15 points an hour 
by the aid of ordinary table machines (Facit ESA and Hamann-Selecta). It was found 
that, except for rounding-off errors in the 8th decimal, even the preliminary solution 


was in agreement with the values in Miller’s table. 


Table I. Start integration for the equation er-— yf by means of successive approximations. The 


computations are arranged as in Table III 1d we will here give only the solution y in the different 


approxin itions in comparison to Miller’s values 


h-y y (2) h-y® 
—0.5 > 898 806 12 898 723 10 10-8 
—0.4 9 541 771 19 541 759 
-—0.3 1) $23 486 55 823 422 39s —.55 823 485 39) 
—0.2 61 787 457 61 787 457 sli 61 787 455 31 
—(.] 4 169 872 67 469 819 25; 67 469 873 25, 
—0 72 901 113 72 901 113 20) 72 901 1138 20 
+0 72 90 113 72 901 113 20) 72 901 113 20 
0 78 106 918 78 106 = 869 16. 78 106 919 16s 
0.2 83 109 270 83 109 265 13 83 109 269 135 
0.3 S7 927 067 87 927 O16 11 87 927 O68 1 le 
0.4 )2 176 =688 G2 976 675 
0.5 7 O72 392 97 O72 337 
There is an error of +1 in the 6th decimal in the first approximation y“ and a corresponding error in the 
Sth de lin the ne pproximation j, The marching process can be based on y@ 
Table II. Calculation of y in the start. The table shows the elimination of the oscillating errors in the first 
approxi! tion by tl n n-value formation, ef. formula (17’’). The difference scheme for F is rather 
irregular in the first approximation, but still it gives almost the same values of 2467 and 2u65F as the 
second approximation. The unity = 10 
First approximation Second approximation 
x Fd oF 2u5F 26h rit oF 26h Q2ueF 
82 163 82 212 
—.3 147 319 147 092 
bo 56 } Zé 64 SU t 222 
2 116 576 7 5381] 116 586 7 337 
51 420 t 260 51 706 3 115 
—.| 93 364 5 376 93 353 5 415 
11 944 1 116 $1 647 2 300 
0 75 528 t O24 75 5385 3 993 
33 584 2 908 33 888 1 693 
61 716 2 955 61 710 2 974 
28 132 17 27 «=822 1 281 
2 50 859 2 260 50 859 2 254 
22 (27 2 213 23 037 973 
3 $2 262 12 262 











YY 
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Table III. A part of the main table for the marching process in the numerical integration of y’ = x — 1. 





x y M+y Y M F —-3VF 








3 —.55 823 4850 5.78 622 388 393 5.78 621 995 —.61 162 615 
995 238 

2 61 787 4584 6.37 266 884 315 6.37 266 569 —.58 176 900 
885 021 

* l -.67 469 8729 6.93 206 008 254 6.93 205 754 —.55 521 838 
792 038 

0 — .72 901 1133 7.46 726 374 206 7.46 726 168 —.53 145 723 

+O -.72 901 1133 —7.11 295 892 —206 —7.11 295 686 —.53 145 723 
—712 939 

. | 78 106 9189 —7.64 066 886 —168 —7.64 066 718 —.51 006 908 
—645 134 

2 83 109 2686 —8.14 735 518 —139 —8.14 735 379 —.49 071 505 
—586 604 

3 87 927 0676 —8.63 500 112 —115 --8.63 499 997 —.47 311 692 
—535 754 

i 92 576 6881 —9.10 532 070 —98S —9.10 531 972 —.45 704 432 
—491 311 

5 97 072 3954 —9.55 980 454 —S86 —9.55 980 368 —.44 230 500 
—452 254 

6 1.01 426 6910 9.99 975 664 —71 —9.99 975 593 —.42 873 737 
—417 755 

7 1.05 650 5902 -10.42 632 411 —61 —10.42 632 350 —.41 620 472 
—387 135 

8 —1.09 753 8455 —10.84 052 100 —53 —10.84 052 047 —.40 459 066 
—359 839 

9 —1.13 745 13815 —11.24 324 799 —46 —11.24 324 753 —.39 379 549 
—335 403 

1.0 1.17 6382 1976 —11.63 530 863 —39 —11.63 530 824 —.38 373 339 
\ comparison with Miller’s table shows that there is an error of +1 in the 8th decimal in the preliminary 


solution; the aftereorrection is so small that it changes only the rounding-off. 


Table IV. Solution of the algebraic equation (26) in the marching process according to the iteration 
£ ] £ £ 


process (28). r gives the order of approximation. 


x=1.0 h=0.1 N = —1.16 353 0863 








r y\? F(r) 

0 — .38 373 343 
] —1.17 632 1977 — .38 373 340 
2 —1.17 632 1976 — .38 373 339 
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Table V. Solution of the algebraic equation (26) in the start by using formula (28) and (31), respectively. 
First approximation. 


xz=0.1 h=090.1 N = —0.76 406 6409 

y Fo) Order y”) Fw 
0 —.53 145 723 
l —.78 178 1650 —.51 118 255 
2 —.78 110 5827 —.51 012 631 ™ —.78 106 8685 —.51 006 829 
3 —.78 107 0619 —.51 007 131 si —.78 106 8685 
"4 — .78 106 8786 —.51 006 845 
5 —.78 106 8691 — .51 006 830 
6 — .78 106 8686 —.51 006 829 
7 —.78 106 8685 


This example shows that the use of (31) can spare a considerable amount of work in the start integration 
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THE TRANSFER FUNCTION OF GENERAL TWO TERMINAL-PAIR 
RC NETWORKS* 


BY 
AARON FIALKOW 
(Polytechnic Institute of Brooklyn) 
AND 
IRVING GERST 
(Control Instrument Co., Brooklyn, N. Y.) 


1. Introduction. Recent years have seen the increased use in many fields of passive 
electrical two terminal-pair networks lacking inductive elements. In particular, those 
networks having a common ground have found widespread employment in servo- 
mechanism design. In these applications, the network function of prime interest is 
generally the transfer function A(p) defined as the ratio of steady state output voltage 
to input voltage in the domain of the complex frequency variable p. 

Our discussion of the transfer function in the sequel will be given in terms of RC- 
networks. However, it is known that by means of simple transformations the results in 
the RC-case can be made applicable to networks containing any two kinds of elements 
only. We will indicate these transformations later for LC-networks, important in classical 
filter theory, and for RL-networks. The actual statement of our results for these two 
cases is left as an exercise for the reader. 

Despite its importance, there have been few investigations of any generality con- 
cerned with the transfer function. Being given a class of networks, the basic questions 
with regard to the transfer function are: (1) by what properties may one characterize 
the transfer functions which arise from this class of networks; (2) given a transfer func- 
tion having these properties, how does one obtain corresponding networks of the class. 


If, as we may, we write 
A(p) = KN/D = K(p" + ap" + --- +4,)/(p” + bp” +--+ + b,) (1.1) 


where .V and D have no common factors, we require a complete characterization of 
K, N and D for the class of networks under discussion as well as a synthesis procedure. 
[t is important to note that since we are interested in what can be obtained from the 
networks alone without amplifiers, transformers or similar devices, the description of the 
multiplicative constant K is just as essential as that of N or D. 

With this view of the problem in mind, we find that the bulk of the literature on 
RC transfer functions, besides concerning itself with quite special networks, only 
partially treats the questions here raised. Thus, generally, the constant K is ignored, 
and incomplete conditions on N and D are given. The complete problem as stated here 
has thus far been investigated in the case of rather specialized classes of networks. 
Bower and Ordung [2] have considered the transfer function of the symmetric lattice 
using a geometric method while we [4] have treated this network analytically as well 
as the L and general ladder networks. 

In the present paper, we do not restrict ourselves to networks of any special internal 

*Presented to the American Mathematical Society, February 24, 1951. (Cf. Bulletin of the American 
Mathematical Society, Vol. 57 (1951) pp. 182-3.) 
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structure but treat the general grounded two terminal-pair (3 external terminals) and 
also the general two terminal-pair (4 external terminals) RC networks (abbreviated 3 
T.N., 4 T.N., respectively.) An informal description of our results follows. For more 
precise statements, see Theorems 1 and 2. It is well known that the transfer function 
of a 3 T.N. or 4 T.N. is a rational function with real coefficients, regular at infinity, and 
having negative real, simple poles. In addition to these properties, we find that for the 
3 T.N., the zeros of N may not be positive real. For given N and D, K must be a number 
in the interval 0 < K < K,, where K, is a constant depending upon N and D which is 
explicitly determined. Conversely, a 3 T.N. may be constructed which realizes any 
transfer function A(p) whose K, N, D satisfy the preceding conditions. In this con- 
nection, Guillemin [5] had previously shown that if a; > 0 in N, then a 3 T.N. exists 
realizing the transfer function A(p) for some unspecified K. In the case of the 4 T.N., 
we find that there is no restriction on the zeros of N. It is necessary and sufficient that 
K lie in an interval —K, < K < K,, where again K, is given explicitly in terms of V 
and D. 

With regard to the above results, we remark that they do not imply the corre- 
sponding complete theorems for any structural sub-classes of 3 T.N. or 4 T.N. such as 
were discussed in [2] and [4]. For the transfer functions associated with a sub-class have 
further distinctive properties which are peculiar to the particular internal network 
structure of that sub-class and which permit synthesis in that structure. 

Finally, we note that the results of the paper may be modified to take account of 
any resistive load or source. 

2. The grounded two terminal-pair network. Theorem 1. The necessary and sufficient 
conditions that a real rational function A(p) given by (1.1) be the transfer function of an 
RC — 3T7.N. are: 

(1) The zeros of D are distinct negative numbers. 
(17) The zeros of N may not be positive real but are otherwise arbitrary. 


(47) m > n. 


(iv) The number K satisfies the inequalities 0 < K < Ko, where Ko is the least of the 
three quantities* K, , b,,/a, , 1ifm = nand of the first two quantitiesifm > n. If Ko # Ka 
then K may equal K, . Here K, is the least positive value of x (if tt exists) for which the 
equation D — «kN = 0 has a positive multiple root. 


Proof: (a) Necessity. As stated in the introduction, conditions (i) and (iii) are well known 
in the case of a 4 T.N. To establish these conditions it is sufficient to write the transfer 
function (1.1) in the equivalent form 

A(p) = ¥Yy2/Y2, (2.1) 


where Y,, and Y.. are the short circuit transfer and driving point admittances, re- 
spectively, of the 4 T.N. and then apply Cauer’s results on RC admittances and Brune’s 
residue conditions. (Cf. [6], pp. 134-136, 211-212, 216-218). Since a 3 T.N. may be con- 
sidered as a 4 T.N. with one input and one output terminal joined together (i) and (iil) 
also apply to a 3 T.N. 

*If a, = 0, omit b,/a, . It may be shown that in any case at least one of the three quantities Kg , 


bm/@n , 1 actually appears as an upper bound for K. 
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lor the proof that the remaining conditions are necessary, the 3 T.N. shown in 
Fig. 1 is considered upon a nodal basis. Here the ground terminal is taken as node 0, 
the “other input and output terminals as nodes 1 and 2 respectively. For our purpose 
the"remaining nodes are identified so that each branch is an Randa C in parallel. Hence 














a ris 


te) 


Fra. 1. 


the admittance y;; (¢ ~ j) of the branch between nodes i and j is of the form ap + 6, 
a = 0, b > 0. Of course y;; = y;; . Write 


t 


4 = _ (@ = 1,2,---, 9d, 


7=0 
i) 


where ¢ + 1 is the total number of nodes. Then y;; is also linear in p with non-negative 
coefficients. Let J be the current impressed on node 1 by the driving source and denote 
by E; the voltage from node 0 to node 7. As is well known [1], the equations of the nodal 


System are 


[= yuk, a Yio, ae eS Yl, 
0 = —Yaky + Yok, — +++ — Yo El, 
(2.2) 
0 = —ynki — Yok, — --> + Yk, . 
Since by definition A(p) = E,/E, , it follows from these equations that 
A(p) = y/o, (2.3) 
where 
You Yos Tee Yor 
— Yai Y33 ee —Y3¢e 
A, = =oap'+eap '+--- +e, (2.4) 
| 
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l 
| Yoo Ys: tee Yee | 
| 7 Ys 8 Yan | 
A, = | = dp’ + dip ‘4 --- +d .d,#0. (2.5) 
Yin Yi: y 


Here A; and A, may have common factors and by (iii) the degree of A, may actually 
be less than s. 
It may be shown that* 


O<e¢<d (j= 0,1, --:,98). (2.6) 


(In order not to interrupt the train of the argument at this point we defer the proof 
to Appendix A). It follows that none of the zeros of A, and hence of A(p) can be positive 
real. This establishes (ii). Also since c; > 0, d, > 0, it follows that in the reduced form 
of the transfer function (1.1), K > 0. Then by (ii), a, > 0 in (1.1). (But note that some 
a, (¢ * n) may be negative). 

In order to establish the necessity of (iv) let us suppose the roles of nodes 0 and 1 
in Figure 1 are interchanged. This new 3 T.N. whose ground terminal is node 1 and 
whose other input and output terminals are nodes 0 and 2 respectively, while all other 
nodes are left unchanged, is termed the complementary network (of the original net work) 
and its transfer function A’(p) is called the complementary transfer function. Then for 
this network under the same impressed current J as before, the input voltage is —£, 
while the output voltage is E, — EF, . Hence 

A'(p) = (E, — E,)/—E£, = 1 — A(p). 


Using (1.1), we find 
A’(p) = (D — KN)/D, (2.7) 
where it is clear that numerator and denominator are relatively prime. Hence the re- 


duced form for A‘(p) corresponding to (1.1) is 
; l | " , 
Pp + ap T oe FF a 
A'(p) = K’ = iF de 2.8) 
a b,p” oop eee + Ob, 
where, as shown above, K’ > 0, ai > 0. Comparison of the two expressions for A’(p) 
given by (2.7) and (2.8) shows that these last inequalities require K < 1 if m = n, 
K < b,,/a, if a, 4 0 respectively. 

The remaining inequality on K in (iv) is obtained by first showing that the transfer 
function of a 3 T.N. has what may be called the “interval property” with respect to its 
multiplicative constant K. That is, if A p) as given by 1.1) zs the transfe r function of a 
3 T.N. then 

A*(p) = K*N/D, 6 < K* < KE, (2.9) 
is also the transfe r function of as TN. 

To prove this, suppose A(p) arises from the network of Figure 1. We may write 
*That d; > 0 is of course well known if either the conductance or capacitance matrix of the y;; 


occurring in Az is positive definite. 
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A(p) as in (2.1). Now modify the 3 T.N. by introducing a new branch joining nodes 0 


and 2 whose admittance is (A — A*)Y../K*. Then the modified 3 T.N. will have as 
short circuit driving point and transfer admittances Y4 = Y..+ (K — K*)Y,../K* and 
Y; Y,» respectively; so that the corresponding transfer function is given by (2.9). 


Now suppose A(p) in (1.1) is the transfer function of a 3 T.N. and let A, be defined* 
as in Theorem 1, (iv). Then K < K,. For by what has just been proved, if K > K,, 
A*(p) = KN Disalsoa3 T.N. transfer function. But then the complementary transfer 
function (D — K,N)/D would have a positive real zero, in contradiction to (ii). 

At first glance the bounds given for K in Theorem 1, (iv) seem quite unrelated. 
However, one may alternatively characterize K, in (iv) as the greatest lower bound of 
all positive « for which the equation D — cN = 0 has a positive real root. It may then 
be shown that this implies the statement as given in (iv). 

b) Sufficiency. Suppose now that A(p) satisfies the conditions of Theorem 1. We 
shall construct a 3 T.N. whose transfer function is A(p). For this purpose, we first 
write A(p) in a special form suggested by the fact that, before algebraic simplification, 
every transfer function has the form (2.3) with d; > c; > 0. We define this special 
form, called an R-function (realizable function), as follows: A real rational function 


R(p (gop + mip ' 4 +++ + g)) (hop' + hp + --- +h,), hy # 0, (2.10) 
is said to be an #-function if the zeros of its denominator are negative real and distinct, 
and if 0 < g; < h; (¢ = 0,1, --- , 1). We eall 1 the degree of the R-function. The ex- 
pression of a rational function as an A-function is not unique. 

[t is shown in Appendix B that a function A(p) satisfying the conditions of Theorem 
1 may be written as an /-function. In the sequel, we give an algorithm for realizing 
every P-function as the transfer function of a 3 T.N. Our method consists of an in- 
duction on the degree / of the R-function. 

The ease | 0 is trivial. If 1 = 1, 1.e., R(p) = (gop + gi)/(hop + h,) with hoh, ¥ 0, 
0< gm <h,,O <q, < h, , then R(p) is realized by an L-network whose series and 
shunt arm admittances are Y, = gop + g, and Y, = (ho — go)p + (hi — g:) respectively. 
Note that the Y,, of this network is hop + h, . 

Now suppose that all R-functions of degree less than / > 2 are realizable as transfer 
functions of 3 T.N. Also assume that for R(p) = U/V of degree f < l, U and V as in 
2.10), the Y., of the corresponding network is of the form V/S where S is a polynomial 
of degree f — 1. We will show that (2.10) is realizable as the transfer function of a 
3 T.N. such that its VY». is of the form >°!_, h;p'*/T where T' is a polynomial of degree 
{— 1 


Let the negatives of the zeros of the denominator in (2.10) bey, <y.2 < ++: <y- 
Choose any o; (¢ = 1, 2, --- ,l — 1) such that 
Ly, Sh KM eS *** See <> 


1 


Then the function Z = ho []i-1 (p + ¥.)/p fe (p + o,) is an RC-impedance who 
canonical expansion is 


Pci +A TS Aete), 4:96. CoOL, ,t<-D Cm 


i=l 


*It follows that Ky is one of the roots of the equation in « obtained by equating the discriminant of 
D — «N = 0 to zero. For the expression of the discriminant of an algebraic equation in terms of its 


coefficients see [3]. 
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In (2.11) the terms of the right member with h, omitted form another RC-impedance 
which we write as c ge (p + &,)/p Pisas (p + ¢;),c > 0. We have 


o <8 < Gy < Se Ste <— *** S Bey KH Cay s (2.12) 


> 


From (2.11) it now follows that 


1-1 — 
h [Tl @t+y) =hp [1] @+o) +e [] @+ 8). (2.13) 
+=] t=] a=1 
If we let 
1-1 t—1 I-1 u ; 
hp I] (p+) = Dd hip’ and cl[pMt+ése = Do hi'p', 
i=1 +=0 i=] ‘=1 
then the h/ and h?’ are all positive; and equating coefficients of like powers of p in (2.13) 
gives 
h=hi, h; = hi + hi’ @=1,2,---,l-—D, h, = hi’. (2.14) 
Since 0 < g; < h,; in (2.10), and in view of (2.14), there exists for each 7 (¢ = 1, 2, ---, 
1 — 1) at least one pair g/ , gt’ such that 
g=gitg, O<g shi, O<g! SH". (2.15) 
Take gi = go and gj’ = g, . 
Now consider the functions 
1 i-1 1-1 I-1 
Rip) = & gp’ ‘/hop T] @ + 0) = D gp’ **/ho TI] @ + @1), 
i i=1 1=0 i=l 
l I-1 
R(p) = 2 gi'p'/e I] (p + &,). 
i=] i=1 
These are R-functions of degree 1 — 1 at most. Hence by the hypothesis of induction 


there exist two 3 T.N. I’, and I, whose transfer functions are R, and R, respectively, 
and whose Y..’s are 


i-1 


1 
YP =h T]@teo)/S, and Y3 =c[[@Mt+eé/.S, 
i=1 t=1 
respectively. Here S, and S, are polynomials of degree / — 2 while \, and X, are arbitrary 
positive impedance-level constants whose values do not affect the transfer functions. 
We shall presently fix \, and A. . 


Now in view of (2.12) we can choose 6; (¢ = 1, 2, --- , l — 1) such that 
&<B, <1 <& << B.<o2 < hhh ms hy 1 <~ Pi S Wes 
Hence 
l 1 i-1 
ys? =hp Tl pt+o)/pt+8), VY? =c [1] @M+t)/Pt+8) (2.16) 
1 1 


t=] 


are both RC-admittances. We are going to modify the networks I’, and I, so that their 


. . . . - /. 1) (2 . 
transfer functions remain unchanged but their Y..’s become Y22 and Y22' respectively. 
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In Fig. 2 consider the networks T'{ and Tj where Z, and Z, are RC-impedances to be 
determined. Evidently the transfer functions of Tj and Tj are the same as those of I, 

















I 2 
om. 2 2’ | e 2 z 
O L, 
/ / 
fi D 
Fie. 2. 
and I’, . However, the Y,.’s of the new networks are 1/(Z, + 1/Y3;’) and 1/(Z, + 1/Y32 


respectively. If these are to equal Y;’ and Y32’ respectively, then it is necessary and 
sufficient that Z, and Z, as given by 
Z, = 1/Yx — 1/Y22, = Z2 = 1/Yn.? — 1/Y22 (2.17) 


be RC-impedances. Now in canonical form 


1-1 i—1 
1/yi2 = B/p+ VB/pt+o), Wye? =Q+ LC/pt+s 
i=l t=! 


with 
B,; > 0, C;>0 @=0,1,---,/—1; 
and 
1-1 pms 
1/¥:? =. DL D./p +o), 1/¥3? =r. Do Ei /(p + &) 
s=1 i=l] 
with 
D, > ©, E; >0 (@= 1,2, ---,t—1). 
Hence 


l—1 t—1 
Z; = B, P + > (B, = \,D;)/(p + o;); Z, a Co + > (C; — 2 ;)/(p + é,), 


and by taking \, and ), sufficiently small and positive, Z, and Z, will be RC-impedances. 
The Y,,’s of the networks Ij, T'; in view of (2.1), are, respectively, 


I-1 Se | 
yi = Y2R = Lo gip'**/ TL @ + 8), 
+=0 t=1 


I 


> g'p'*/ I] (p + B,). 


i=1 


(2 (2) 
Yio — Yoo R, 


Now connect the networks T{ and I} in parallel to form a new 3 T.N., [. Then 
(1) 


[6 p. 146] the Y,. and Y,.,. of T are given by Yi2 = Yio + Vis, Yo = Ys + Ys 
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Replacing the Y,.’s and Y2.’s by their expressions given above and using (2.13) and 
is (2.10) while 


Y 


(2.15) we find that the transfer function A(p) of I 


I-1 


too = > hip I] @ + £6;). 


i=1 


This completes the induction*. 

The preceding synthesis procedure may be modified to accommodate any given re- 
sistive source and load. In the general case, this may require a finite reduction in the 
value of the maximum realizable K. However, if the source is zero or the load is infinite, 
all K except possibly AK, may be realized. 

As an illustration, we outline the synthesis procedure for the case of a resistive load, 
which we may take as 1 ohm without loss of generality. If A < K, we may suppose 
the R-function in (2.10) which corresponds to the transfer function is such that 0 < 
g; < h. . The procedure is now started by choosing 7) ae l) such that 0 < 
% KX % s Te < *** rT, < y,. Then as before we have 

l 


I I . 
h T] ee oF =h I] @+ r) +e [1 @+x, ie | D 


t=1 


where 7; < xi < Te < x2 <*** < x:-1 < 7. The original synthesis procedure is now 


applied to 


to give a 3 T.N. T whose Y wc is 


i 1 
ull @t T;) [l@+x), p > 0. 

If the 7; are taken sufficiently close to the y; this will always be possible. The network 
consisting of I working into the load then realizes the original transfer function. 

We conclude this section with the observation that theorems analogous to Theorem 1 
exist: for 3 T.N. having only resistance and (self) inductance er capacitance and (self) 
inductance. As shown in [4, Sec. 7] these are obtained from Theorem 1 by replacing p 
by 1/p in the RL case and by replacing p by p* in the LC case. 

3. An illustrative example. As an illustration of the foregoing theory and s 
procedure consider the function A(p) = KN/D = K(p° — p+ 9)/(p + 1)(p + 14). 
First let us determine for what range of K, A(p) is the transfer function of a 3 T.N. 
Since m = n = 2, a, = 9, b, = 14, Theorem 1 tells us that 0 < K < Ky where Ky = 
1]. To ealeulate K, it is best to proceed indirectly by first eliminating 


ynthesis 


min [K, , 14/9, 
x between the equations D — «N = 0 and d/dp (D — «N) = 0. This gives 16 p> + 


10 p — 149 = 0 whose roots p, = —3.380, p. = 2.755 are the multiple roots of D — 
«kN = 0 which occur for all variations of x. Only the latter of these is positive and corre- 
sponding to it we find «x = D(p.)/N(p.) = 4.546 = K, . Hence Ay = min 


(4.546, 14/9, 1] = 1. 


*A variation of the above procedure may be obtained by modifying the decomposition as given in 
(2.13). Using either procedure an appropriate choice of (2.15) will frequently lead to simple network 
realizations. Still another modification of the synthesis procedure makes it possible to always choose 


Z; as a condenser and Z, as a resistor in Fig. 2 at each stage of the synthesis. 
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Let us choose (as we may by Theorem 1) K = K, = 1. To apply the synthesis 
procedure, we must first write A(p) as an R-function. Using the ideas of Appendix B 
(or by inspection) we find that the introduction of the common factor (P + 3) into 


A(p) converts it into the R-function 

R = (p° + 2p? + 6p + 27)/(p + 1) (p + 3) (p + 14). 
Following the procedure of §2(b) choose o, = 2, ¢. = 6. Then (2.13) reads p* + me! + 
59p + 42 = pip + 8pt+ 12) + (10p° + 47p + i: from which £, = 6/5, & = 7/: 
Now, for simplicity, in (2.15) let 2 = 2 + 0, 6 = 0 + 6; also choose 6B, = 3/2, 
Then we get 


R, = (p* + 2p’)/p(p + 2)(p + 6) = p/(p+ 6), R. = (Gp + 27)/10(p + 6/5)(p + 7/2), 


2. 
B. = 4. 


while according to (2.16) 
Ys. = p(p + 2)(p + 6)/(p + 3/2)(p + 4) 
and 
Yo. = 10(p + 6/5)(p + 7/2)/(p + 3/2)(p + 4) 
respectively. 
The transfer function R, being of first degree is realized by an L-network IT, whose 
series and shunt arm impedances are Z,, = A,/p and Z,, = ,/6 respectively, and whose 


Y3, = (p + 6)/A, . Use (2.17) to determine Z, taking \, = 3/8 for simplicity. Then 
Z, = 1/2 p + 1/8(p + 2). Finally form Ij as in Fig. 2 


\s for R, , it is of the second degree and hence the reduction procedure is repeated 
for it to obtain R, and R, of the first degree. We give the results, the notation being 
evident from what has preceded. 


Op + 27 71 g, By o; Y2 


hs = 10 + 6/5\(p + 7/2)’ 6/5, 14/9, 7/4, 2, 7/2” 
l0p° + 47p + 42 = 10p(p + 2) + (27p + 42); 
Rp. =0 yi. 1Op(p + 2) . R, = 6p + 27 | itt se 27(p + 14/9) | 
oe ‘ — (p+ 7/4)’ ‘4 27(p + 14/9) ’ ve (p+ 7/4) 


The transfer functions R; and R, are realized by L-networks I , I’, whose series and 


shunt arm impedances are Z;, = ©, Z;, = (p + 7/4)/10p(p + 2); Z,. = Ax/(6p + 27), 
Zs, = /(21p + 15) respectively; and whose Y..’s are Y3; = Ys* and — = 








Fia. 3. 
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27(p + 14/9)/d, respectively. Hence we need only determine a Z, using (2.17) and taking 


\, = 7/36. This gives Z, = 1/27. Now form Ij as in Fig. 2 and denote by I, the parallel 
combination of T, and Tj . Introducing the impedance level factor A, in [, the Yo. of 
I, is given by ¥2.’ = 10(p + 6/5)(p + 7/2)/d2(p + 7/4). Again use (2.17) to determine 
Z, taking \. = 4/7. Then Z, = 1/10 + 4/175(p + 6/5). Forming T; from [, and Z, 


as in Fig. 2, the required final network is the parallel combination of '{ and T'; as shown 
in Fig. 3. 

4, The two terminal-pair network. Theorem 2. The necessary and sufficient conditions 
that a real rational function A(p) given by (1.1) be the transfer function of an RC — 4 T.N. 


are: 


(i) The zeros of D are distinct negative numbers. 


(22 ma nN. 

(772) The number K satisfie s the mnequalitre s —K, < K < K,, where K,j is the least of 
the three quantities Rat i sent. 9 if m = n and of the first two quantities ifm> n. 
If K, # | Ky) then K may actually equal +K, . Here K, is that real value of « of smallest 
absolute value (if it exists) for which the equation D — kN = 0 has a positive multiple root. 


Proof (a) Necessity. For (i) and (ii) see the remarks in Sec. 2. To prove (iii) consider 


the 4 T.N. of Fig. 4 on a nodal basis taking the input terminals as nodes 0 and 1, the 














A e 
a" = 
oo 3 











Fia. 4. 


output terminals as nodes 2 and 3, and choosing the remaining nodes as in the 3 T.N. 
case. Let E,, and FE, be the input and output voltages, respectively. Then in the nota- 





tion of Sec. 2 we have A(p) = Ey./E,, = (E, — E;)/E, . It follows from (2.2) that 
A(p) = (A, — As;)/A: (4.1) 
with A, and A, as in (2.4), (2.5) and 
Y31 Y32 Y34 Yse | 
| — Yaoi Yoo — Yo vie —Yae | 
| . ,_.8- 
A; = | Ya — Ye Ya le —Ysr | = Cop + C\p +e + 
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Since A, is of the same form as A, we have 


0<c<d (j=0,1,-:-,9, (4.2) 


in addition to (2.6). 
If we equate the two expressions for 1 + A(p) obtained from (1.1) and (4.1), we 
have 


(D + KN)/D = [A, + (A, — As)]/Ae « 


Since the fraction on the left is evidently in reduced form, and in view of (2.6) and (4.2), 
it follows that | K | < 1 if m = n and that D + KN = O has no positive real roots. 
This in turn implies that | K | < | b,,/a, | if a, # 0, and that | A | < | K,| for Ky as 
defined in Theorem 2. This last fact follows as in the 3 T.N. case; for the “interval 
property” of K holds fora 4 T.N. with 0 < | K*| < | K | in (2.9). 

(b) Sufficiency. Let A(p) satisfy the conditions of Theorem 2. Then proceeding as 
in Appendix B we find that both D + KN and D — KN have leading coefficients 
positive and have no positive real zeros. Thus there exists a polynomial P having distinct 
negative real zeros such that P(D + KN) and P(D — KN) have non-negative co- 
efficients. We may suppose P and D relatively prime. If 


i l 
PD= > hp'' and KPN= > gp’, 

i=0 t=0 
it therefore follows that |g,;| < h; @ = 0,1, --: , D. Write A(p) = KPN PD = 
(N, — N.)/PD, where N, consists of those terms of APN having positive coefficients 
while —N, consists of those terms of KPN having negative coefficients. Then V,/PD 
and N./PD are R-functions. Hence by the results of Sec. 2 there exist two 3 T.N., 
r, and I’, whose ground, input and output terminals are 0, 1, 2 and 0, 1, 2’ respectively, 
and whose transfer functions are A, = N,/PD and A, = N./PD respectively. Now 
form the 4 T.N. shown in Fig. 5 whose input terminals are 0 and 1 and whose output 


— 





























| 
— | 02! 
2 
0 
Fia. 5. 
terminals are 2 and 2’. The transfer function of this network is evidently A, — A, = A, 


temarks similar to those of Sec. 2 may also be made here concerning source and load, 


and the RL and LC networks. 
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APPENDIX A 


Consider the determinant 








Ras Aye oe eee Bin 
— 21 A22 — oz em — Gon 
| 
H = | —@31 — A32 433 ike — A3n |, 
— An — Ano —An3 see a. 
where a,, = >.” z; 4;;(¢ = 2,3,---,n) and all the a,;(i ¥ 7) and a,, are independent 


pe 
variables. We shall prove that H is a polynomial (multilinear form) in the a;;(t- # Jj) 


and a,, with non-negative coe fhicre nts. 





This is evidently true for n = 1, 2. Suppose it is true for determinants of the form 
H of order less than n(n > 3). We may write H = H, + H,.where H, contains those 
terms which are independent of a,,(i = 1, 2, --- , m — 1) and H, those terms which 
are linear in each a;,(i = 1, 2, «++ , nm — 1) separately. In H, if we let a;, = 
O(c = 1, 2, --- ,n — 1) we find 
| Ay; a, Qj (n-1) 
| | 
| —doy as — 2 (n-1) 
i. = Ann | ’ 
| saad bee aceh in cunkeaken wrneae eT 
| , 
—Q(n-1)1 — GQ (n-1)2 Q(n-1) (n-1 
where a’, = a,; — a;,(i = 2,3, --+ ,n — 1). Hence applying the hypothesis of induction 


we conclude that the terms of H, have positive coefficients. 
As for H, each of its terms is contained in at least one of the polynomials a,,, dH /da,, 
(= 1,2, ---,n— 1). Now 





| . n 
| O 0 0 cee l 
| 
. | —o Ao2 —Ag3 see —Qon 
oH | zi 22 74 
Oa, a 
NR cl alla a al ttt thal aati 
on™ —An2 —da Qny 
| an a, Ons Bn(n-1) 
i —@ Ao —A>2 “ — Qo(n-1) 
Q(n—1) (n-1 
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and fori = 2,3, ---,n-—1, 

| Qi Aye ar ai; — Qin 

[ =e a — Ao; sl als on 

| 

OH 
O0;n - 
0 0 l see ey | 

| TrTUTCTEUCU CTCL TTT ee ee eoeveccece 

| —~—Aa1 —An2 eran —Ani — Onn 

ay Aye siecle Q4(i-1) Qy (441) sl (a; + 4,,) 

— 2) Az2 28° —A2(i-1) — A2(4+1) -++ —(@2; + G2.) 

= — ai, 1s ~"Gu-12 °°" Q(i-1) (4-1) ~@i-1) (541) ° °° —(a;, 16H Ac5-1)0) . 

| —Gc¢41)) ~Auene °° * O40 G-2) Qii+iy(ita) °°* —(Qci41ys HH Aci siya) 

| 

| | 

—~— O41 — Gao sic —~ Bas 1) =~ On(i+1) sata (Aan — a,;) 

Replacing a,) + a;,(¢ = 2, 3, , n — 1) by new variables a{, in the first case, and 
a;, + a;,(j7 = 1, 2, --- , — 1) by new variables a{, in the second case, we see that 
the final determinant in each case is reduced to a determinant of the same form as H 
but of order n — 1. Hence the terms of H, also have positive coefficients and the above 


italicized statement about H is proved. 

Now make the following specialization: n = t — 1, @4. = Yasin = 1,2, ---,é— 1), 
Q;;5 = Yuanugen(t = 1,2, «++ ,t — 17 = 2,3, --+ ,t — 1 t F J), Gio = Yura + 
Yorsryo(@ = 1, 2, +++ , £ — 1). Then H becomes A, and we find that A, is a polynomial 
in the y,;(¢ ¥ j) with non-negative coefficients. This proves the c; > 0(7 = 0, 1, --- , 8) 
in (2.4). 

In A, replace the elements of the first column by yoo, — Yao, *** » — Yeo respectively, 
to form Af. Then by symmetry A{ is also a polynomial in the y,;(¢ # Jj) with non- 
negative coefficients. Now in A, of (2.5) add the sum of columns 2, 3, --- , 4 — 1 to 
column 1. We see that A, = A, + Aj. This proves that d; > c;(j = 0,1, +++ , 8). 


APPENDIX B 


Lemma*: Let F(p) = ap* + --+ ,a@ > 0 beareal polynomial having no positive real zeros, 
Then there exists a polynomial P = [](p + 4,) with the 6; > 0 and distinct, such that 
PF is a polynomial with non-negative coefficients. 

*For a more general result which includes this one see [7]. We include the following simple con- 
structive proof so as to keep the paper self-contained with regard to the synthesis procedure. 
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The real irreducible factors of F are either of the form p+ +7 with r >0 or 
p — 2p cos gp + p with p > 0,0 < ¢ < 2z. It clearly suffices to prove the Jemma for 
F equal to this latter quadratic polynomial with 0 < cos g < 1. Replacing p by pp we 
q 1 pol} I P DY PI 


may without loss of generality suppose that p = 1. 
Let P; = (p + 1)’ where 0 is an integer such that 
b > 2 cos ¢g/(1 — cos ¢g). (B.1) 
We have 


Pe = > 6p’, 


where 


/ .\ 
a.= (") —2 cose. ° ) 4. (_ b J (4 = 0. i arn , 6b 4- 2), 


fr \ 
f \ 


and | .) is the binomal coefficient, with ) defined as zero if 7 > bor7z < 0. 


Now 0 = 6:42 = 1> 0; and 0, = 4,.,; = b — 2cos¢g > Oby (B.1). If 2<i< 3B, 
b (b—171+1)(b —i1+2 (4 — 
9, = ( ) : J - U pee a ere | (B.2) 
1 - e ’) b — 7 + 2) 
We have 7(b —74+ 2) < (¢+6—7-4 2)°/4 = (b+ 2)°/4;and (b —7+1)(b-—7+2) + 
4 — 1) = {b(b + 2) + [b — 27 — 1)]°}/2 > b(b + 2)/2. Hence in (B.2), the bracket 
satisfies the inequality, | ] > 2b/(b + 2) — 2 cos g > O by (B.1). Since the co- 


efficients 6; of P,F are all positive we may by continuity considerations form P{ = 
[Ti. p + 6;) with the 6; sufficiently close to 1 and distinct, such that P{F has positive 
coefficients. This proves the lemma. We note that the method used in the above proof 
may not necessarily lead to the simplest polynomial P. 

Now let A(p) satisfy the conditions of Theorem 1. Then evidently KN is a poly- 


nomial such as F in the lemma. But D — KN is also such a polynomial! To prove this 
we first show that D — KN has no positive real zeros. Consider the zeros ¢;(x) of D — 
xN = 0 where « is a non-negative, real parameter. When x = 0 these zeros are the zeros 


of D and hence are negative real and distinct. Then for small positive x the ¢, will be 
negative real and distinct. As x increases positively, positive real ¢; may arise only if 
one of the following has occurred: (a) a ¢; has gone through zero, (b) a ¢; has gone 
through ~, (c) two or more complex ¢; have become equal on the positive real axis. 
In case (a) we must have x > b,,/a,(a, # 0); in case (b) «x > 1 if m = n; in case (c) 
x > K,. Hence, by (iv) in Theorem 1, if « = K, D — KN has no positive real zeros. 
We must still show that the leading coefficient of D — KN is positive. This is evident 


unless m = n, K = K, = 1. In this case, we have D — N = p >a (b,-i — Gn-i)p’, 
bn» — In-» * 0, r < m. Suppose the leading coefficient b,,_, — a,-, < 0. Consider 
the polynomial Q, = tte 0 €m-iP', Where ¢€,-; = b,-; — Gn-; if a,-; < 0 and e,-_; = 
b,-; — (1 — e)an-:, > Oif a,_; > 0. For e, sufficiently small the leading coefficient 


of Q; will be negative so that there exists a p) > 0 such that Q,(po) < 0. Now 
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D(po) — (1 — &N(po) = «€ >. An-iPo + »» [bn-«+ — (1 — ea,-,]po 
t=r+l 


t=0 


<< } An—iPo + Qi(po) for <<a. 


t=r+l 


Hence for « sufficiently small, say « = e, < €, , D(po) — (1 — &)N(po) < 0 while 
D(o) — (1 — &)N(@) > 0. Thus D(p) — (1 — &)N(p) = 0 has a positive real root, 
which contradicts what has been proved above about D — «kN = 0 for « < K,. Thus 
by the lemma, there exist two polynomials P, and P, such that KP,N and P,(D — KN) 
have non-negative coefficients. Then P,P, will simultaneously convert KN and D — KN 
into polynomials with non-negative coefficients. We may suppose P, , P, and D are 


relatively prime in pairs. Now let 


I l 
A(p) = KN/D = KP,P.N/P,P.D = > gp'‘/ Do ho’. 
+=0 i=0 
Then the zeros of the denominator of the last fraction are negative real and distinct, and 
by considering KP,P,N and P,P,(D — KN) it follows that 0 < g; < h;(é =0,1,---, D. 
Hence the last fraction is an R-function. 

We remark that if Ky = 1 or b,,/a, and if R, is an R-function corresponding to Ky, , 
then AR,/K, corresponds to any K in the range 0 < K < K, . However, if Ky = Ky, 
no such R-function Rp, exists and the degrees of the R-functions corresponding to K in 
the range 0 < AK < K, must increase indefinitely as K — K, . 
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A NON-PLANAR BOUNDARY PROBLEM FOR THE WAVE EQUATION* 


BY 
GEORGE K. MORIKAWA 


Hughes Aircraft Company, Research and Development Laboratories) 


1. Introduction. Notwithstanding the relatively recent outburst of activity on the 
part of aerodynamicists concerning problems dealing with the wave-equation, some 
problems with radiation type conditions still deserve attention. An important class of 
problems of this type are those with non-planar, non-axially-symmetric boundaries, in 
particular, wing-body problems in linearized supersonic flow. For slender wing-bodies 
[1], where variations in the flow direction are small, the stationary three-dimensional 
problem is reduced to a two-dimensional potential (Laplace) problem in the cross-flow 
plane. For conical wing-bodies [2], where the surface is generated by linear elements 
passing through the vertex, a reduction from three to two-dimensions is evident, and a 
ransformation (Chaplygin) further reduces the subsequent elliptic problem to a two- 
dimensional potential problem. However, for more general type boundaries, the problem 
is relatively difficult. One of the first attempts in this direction has been made by 
Ferrari [3], using an iteration procedure which is also discussed by Lagerstrom and 
Van Dyke [4]. In this paper the problem, which the author considers to be the funda- 
mental wing-body problem for stationary linearized supersonic flow, is discussed. An 
approximate solution is obtained and expressed in terms of the pressure, a quantity of 
interest to aerodynamicists. 

2. Formulation of the problem. Consider the wing-body configuration shown in Figure 
1, where the body is an infinitely long, circular cylinder (normalized to radius =1) 
and the wing is a semi-infinite flat plate with the leading edge normal to the free-stream 
direction defined by the velocity W, i.e., zero sweep-back. W is the uniform flow (nomi- 
nally in the z-direction) far upstream and away from the body and the dashed lines 
from the leading edge-body junction are elements of Mach cones indicating supersonic 
leading edges. The problem considered here is the case where both body and wing are 
at a small angle of attack, a, with respect to the uniform flow at — ©; the more general 
case, where the body is at a and the wing has an arbitrary angle of attack distribution, 
can be handled as easily. Since the leading edges of the wing are supersonic, the flow 
field upstream is known everywhere and the boundary conditions are of the first kind [5], 
where one component of the velocity, namely, v = g, , is known everywhere in the 
plane of the wing (2,z plane) and on the body. By a simple super-position procedure 
made possible by the linearization, the original problem with both body and wing at @ 

*Received July 19, 1951. This paper is part of a thesis submitted in partial fulfillment of the require- 
of Doctor of Philosophy at California Institute of Technology, June, 1949, and a 
subsequent report, “The Wing-Body Problem for Linearized Supersonic Flow”’, Jet Propulsion Labora- 
tory, California Institute of Technology, Progress Report No. 4-116 (Dec. 19, 1949), written under 
U.S. Army Ordnance Contract No. W-04-200-ORD-455. A more complete description of the physical 
problem (and related aspects) than given in this paper is presented in the report, hereafter referred to as 
the “JPL report’. The author wishes to thank Professors H. J. Stewart, P. A. Lagerstrom, A. Erdelyi 


and C. R. DePrima for their suggestions and criticisms. 
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can be replaced by an equivalent problem* with body at zero angle and wing at angle of 
attack a, plus an apparent twist imposed by the original flow upstream of the leading 


vr 
x@—— MACH LINE 








"oO 











Fic. 1. Wing-body Configuration. 


edge. Then the boundary conditions in cylindrical polar coordinates, in terms of the 
perturbation velocity potential, are: 


g,(7,0,z2) = — Wal + 7} > %; y-(1,0,z) = 0, (la) 
, 


where ¢ satisfies the wave equation 


l 1 
Ay —- 0., = yg + : Y; a r 09 —~ Goze = 0. (1b) 
/ 
For arbitrary Mach number M, the coefficient of Yrs is —(M’* — 1), but there is no 
loss in generality in taking M = 2’, since, by a simple similarity transformation [5] 


(linear), the arbitrary M case results. The initial condition (z is a time-like direction) 


for the equivalent problem is taken to be 
g(r,6,0+) = ¢.(r,0,0+) = 0. (le) 


Actually, ¢.(r, 0, 0+) ¥ 0, that is, g, # 0 if z = O is approached from positive z in 
the x, z plane, but this condition causes no difficulty here.** Then this radiation type 
problem given by Equations (la), (1b) and (1c) is a natural problem for the methods 
of the Laplace transformation;T and subsequently it will be evident that these methods 
lead to an approximate solution of an inherently complicated problem. 

*Actually, the equivalent problem consists of two problems, the second consisting of the body alone 
at angle of attack a; but since the pressure component ¢, is the quantity of usual interest, the second 
problem contributes nothing. 

**A necessary condition for uniqueness is the condition of outgoing waves. 

7Several British workers have recently applied Laplace transform methods to supersonic flow 

problems; see, e.g. [6,7]. 
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Assuming that the solution is twice differentiable and that the second derivatives 
have Laplace transforms, the equivalent problem given by Eqs. (1) becomes the trans- 
formed problem 


Ay-—sy=0 (2a) 
with the boundary conditions 
] Wa 1 
~ Yolr,038) = ——*(1+5), — v,(1,8) = 0, (2b) 
7 s F 
where 
Vir,6;8) = | e“elr,0,2) dz (2c) 
“0 
is the Laplace transform of ¢ with respect to the free-stream direction, z, and s is the 
transform variable. (A convenient notation is to write Y = L{¢; s} and the inverse 
Laplace transform ¢ = £ '{y; z}; then the pressure p = —plWg, , where p is the free- 
stream density and g, = L ‘{sy; z}). For convenience in presentation, since none of 


the essential features are lost, the sweepback is taken equal to zero. The differential 
equation (Eq. 2a) is the two-dimensional modified Helmholtz equation with the param- 
eter s(Re s > 0) and the transformed problem is an elliptic boundary-value problem in 
the r, @ plane. For strong conditions [8, 9] on y exterior and on the boundary of the region, 
i.e., W is twice continuously differentiable in (r, @), the solution is expressible as an integral 
in terms of ¥ and the outward normal derivative y, on the boundary by application of 
the Green’s theorems. But y and y, are related on the boundary C and the solution 


takes the form 


ures) = -2[ (v%-o%)a, (3a) 


7 
2r on on 
where G is the sum of a singular and a regular function in the region. If G is chosen 


such that dG/dn OQon C 


1 0 . 

ieen =u +f oa (3b) 
2r Jc (On 

where G is called the second Green’s function. The Green’s function is then characterized 

by the following conditions: 


1) G(r,@;r’,0’) = G(P;Q) satisfies Equation (2a) for P(r,é) ¥ Q(r’,6’); 


2) G(r,@;r’,@’) has the proper singularity at P(r,@) = Q(r’,@’); (4) 
IG 

3) : : = QonC. 
on 


A consequence of the first two conditions is that G(r, 6; r’, 6’) = G(r’, 6’; r, 0) symmetric 
in (P; Q). Weaker conditions, e.g., the condition that y, may be discontinuous, but 
integrable on the boundary, are permissible for Eq. (3b) as in potential theory [10]. 
The fundamental solution of the modified Helmholtz equation (Eq. 3b), i.e., the 
solution independent of @ and singular at the origin, is K,(sr), the modified Bessel func- 
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tion of the second kind of zero order [11]. For a fixed point P(r, @) with a variable point 
Q(r’, 6’), the fundamental solution may be written immediately as Ay(sp) since the 
differential equation is invariant under translation; p = [r? + r’” — 2rr’ cos (@ — 6’)]'” 
is the distance between P and Q. This solution may be interpreted as the Green’s function 
for the entire z, y plane, since it satisfies the required properties (Eq. 4). 

In determining the Green’s function for a region with a given boundary, the function 
sought will always have the form 


G(P; Q) = Ko(sp) + A(P; Q), (5) 


where /(p; Q) is regular in the region. Since Ky(sp) is symmetric in (P; Q) then h(P; Q) 
must also be. The invariance properties of the differential equation will be helpful in 
determining h(P; Q) and give some intuitive meaning to it. The Helmholtz equation 
(Eq. 2a) can easily be verified to be invariant under the following transformations: 


1) translation, (7, 7) = (4 + a,y + 5), 
where a and b are constants; 

2) reflection on the axes, e.g., on the x-axis, (6) 
2, ¥) = (zx, —y) or reflection on any straight line; 

3) rotation, (7, @) r,@-+ c) where c is a constant. 


There appears to be no simple invariant transformation with respect to inversion on the 
unit circle, as there is for Laplace’s equation. Such a transformation would be helpful 
in obtaining an intuitive notion on constructing the Green’s function for the circle [12]. 
For the transformed wing-body problem given by Eqs. (2) the solution can be con- 
fined to the upper half-plane and the complete solution obtained later by using the 
known symmetry conditions on the solution with respect to the x, z plane. Then the 
boundary is that given in Fig. 2. The Green’s function for the upper half-plane is 


G,(P;Q) = Ko(sp) + Ko(sp,), (7a) 


where 


p, = [r? +7” — 2rr’ cos (6 + 0’)]'”” 


is the distance between P and the reflection of Q on the z-axis. The Green’s function 


for the unit circle is easily obtained by using the addition theorem [11] for Ay(sp): 


K,(sp) = I(sr)Ko(sr’) + 2 2, [,(sr)K,,(sr’) cos [n(@ — 6’)], rar (7b) 
with r and r’ interchanged for r > r’, where J, and K,, are modified Bessel functions 
of the first and second kinds, respectively, of the nth order. 

For r < r’ consider h(P; Q) of the form 


h(P;Q) = A,K(sr) + >. A,K,(sr) cos [n(@ — 6’)], (7¢) 


where A, and A, are determined to satisfy the boundary condition 
aG.(P;Q) | 


or ir=1 


= 0; 
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then 


_ Ii(s) 
Ki(s) 


I*(s) 
=a I fey” 7d 
Ki(s) )s (7d) 
where the primes on the Bessel functions denote differentiation with respect to r and 


(8) ge / 
G.(P;3Q) = K,(esr’ E (sr) — Los Ke | 
L K,(s) 


A, = K,(sr’), A, = —2 


’(s 
+2> Ko] 109 — a : Kc | cos [n(@ — 6’)], yr <r (7e) 


with 7 and 7’ interchanged for r > 7’. 


P(r,@) 





CONTOUR C 
TITTIT TTI Ea Es EERE ES ml 
‘ 
Q(r',-6) 





Fic. 2. Boundary for Transformed Problem. 


The Green’s function for the boundary given in Fig. 2 is obtained simply by using 
the reflection property on the x-axis with the results for the circle (Eq. 7e) and 


;, (8) 
G P;Q) = IK wr) To(or — To(8) Keo | 
K,(s) 


i. 


+4) Kw] 1467 — an Ks) cos n@ cos né’, Peer (8a) 


with r and r’ interchanged for r > r’. A convenient equivalent form to Eq. (8a) is 


, . . is) .. . 
G(P;Q) = K,(sp) + Ko(sp,) — of se K,(sr) K,(sr’) 
K,(s) 


> 7 
+ 2 pe I,(8) K,,(sr)K,,(sr’) cos n@ cos no’ | (8b) 
n= K,(8) 

With the determination of the Green’s function the transformed problem is formally 
solved. There follows, for each solution of specific problems, the interpretation by the 
inverse Laplace transformation into the solution of the original physical problem; some- 
times this is a formidable task. 

Some intuitive meaning may be given to the Green’s functions for the transformed 
problem. For example, the inverse Laplace transform of the fundamental solution 
K,(sp) represents a supersonic source singularity in the unobstructed physical space. 
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3. Splitting of the transformed solution y. Since on the circle y, = —y 0, the 
transformed solution by Eqs. (2b) and (3b) is 


Wa f° 1 
wir. G38) = = = | (1 + 53 )(G(r, 837’ x) + G(r, 0;r’,0)| dr’. (9a) 
1 


27s ? 


Using the second form of the Green’s function (Eq. 8b) 


Wa ff [ 1 7 
¥(7r,0;s) = - | 1+ LY Kols( + r’° — 2rr’ cos 6) "| 
Jy yy 


TS 


1 Poe “ p2 ¢ , , 21 , 2Wa 4 ' l [‘(s . P 
- K,[s(r -— fT -+ 2rr’ cos 6’) 1} dr’ — - | l — “35 =~) K, (sr) Ko (sr ) 
Ji rt Bold 


Ts 


+2 =. ee K,,,(sr)K2,(sr’) cos 2n6 | dr’ 


yo+y”. 9b) 


This transformed solution is confined to the first quadrant; the solution in the other 
quadrants are obtained by symmetry arguments. The solution given by Eq. (9b) is 
written as two integrals, y'"’ and y*~’, respectively, for two reasons: 1) each can be 
given a physical interpretation, and 2) the difficult core of the solution, i.e., the second 
integral y'*’, is separated from y’’, the inverse Laplace transform of which can be 
obtained simply. By inspection, it is evident that y ” can be regarded as the transformed 
solution for the equivalent problem (cf. first paragraph of Sec. 2) where the circular 
cylinder (a = 0) is replaced by a flat plate (a = 0); thus y"” is called the “flat plate” 
solution. Then the second part of the solution y~’ can be interpreted as the correction 
solution needed to satisfy the boundary condition of zero flow through the circular 
evlinder, i.e., the circular cylinder is a stream tube; thus y*~’ is called the “body”’ solution. 
The principal difficulties arise in attempting to perform the inverse Laplace trans- 
formation of yy” 

4. The ‘‘flat plate” solution ¢"’’. The “flat plate” solution in terms of the pressure 
component (cf. foot note following Eq. 2c) ¢.'’ on the surface is conveniently written 
in Cartesian Coordinates. The transformed solution for gS"? = 2° '{sp""’; z} by Eq. (9b) is 


Wa f oe ne 2)1/24 
sy f.u°s) = pa | 1+ 2) {K as — € = a y) } 


a“ j \ bad 


+ Kyls((2 + 8? + y*)'7]} dé. (Ge) 


The order of integration can be interchanged and the inverse Laplace transformation 
is performed first: 


£'{K,[s((z — &)? + y’)'”7]; 2} 


= 0, tE< [x — (2 — y’)'”] and [e+ (2? —y)'?] <é (10a) 
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= 0, (22 — y)'*® —a] <é (10b) 


= [2° _ y sa + t)”] _ 0 <ét < [(2 = y’)'” _ x] 


In addition, since the lower limit of the integral is 1, three regions of integration are 








MACH LINE 





Fria. 3. Regions of Influence in the z, z Plane. 


defined and have a physical interpretation (ef. Fig. 3, where these regions are shown 


in the x, z plane): 


Region I, (x — z) > 1: Region of influence of the Mach cones from the wing leading 
edge 
Region Il, —1 < (x — z) < 1: Region of influence of the Mach cone from the wing 


leading edge-body junction plus Region I. 


fegion III, (2 — z) < —1: Region of influence of the Mach cones from the opposite 
wing leading edge plus Region I plus Region IT. 


The ‘flat plate’ solution, ¢:'’, on the wing surface (z, z plane) in these three regions is 
given in terms of elementary funetions (cf. JPL report) and is shown graphically in 
Fig. 4 by dashed curves. For the chordwise pressure distribution (Fig. 4b) attention is 
focused primarily on the body (a fairly complete description of the solutions is given in 
the JPL report). 

An additional solution which can be obtained immediately is important in the dis- 
cussion of the complete solution. The “flat plate” configuration (cf. Sec. 3) is modified 
now by inserting a semi-infinite plane barrier at x = 1, parallel to the y, z plane down- 
stream from the leading edge-body junction. This modification does not affect the 
solution g‘'’ outside the Mach cone from the leading edge-body junction; but inside 
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body with infinite radius. Thus, in this region, the solution (called the “modified flat 






this region of influence for x > 1, the inserted plane simulates the limiting case of a 
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Fig. 4. Pressure Distribution on Surface: A) Spanwise, B) Chordwise. 


plate”’ solution g* 
complete solution. The pressure component is given by 


; a 
“1 \ s 


g*(z,0,z) = . [ ~~ (; + te — (x — §)*]"'”’ dé 


) can be interpreted as an upper bound (at least in Region II) to the 
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and is given in Fig. 4 by the light solid lines. Note that the first integral is the “flat 
plate” solution in Region II (cf. JPL report) and the second integral can be interpreted 
as the correction solution needed to satisfy the boundary conditions of zero flow through 
the semi-infinite plane barrier (ef. interpretation of “body” solution in Sec. 3). 

5. Remarks on the “‘body” solution y*’. From the comparison of the “flat plate”’ and 
‘“‘modified flat plate’’ solutions in the previous section and the discussion on the splitting 
of the complete solution in Section 3, it is evident that the “body” solution, in terms 
of the pressure component ¢,”’, is effective only within the Mach cones from the wing 
leading edge-body junction and is identically zero outside this region. The inverse 


Laplace transform required in the particular case of pressure on the body, r = 1, is 
(K.. (er! Te 
| Ko,(sr’) 1 [ .| K.,(sr’) 
2 Na 2; = —— | ¢ 1; Ide rf? Lest 12a 
| sK;,,(s) \ 2mi Jy SKoz,,(s) . ps acs iit (12a) 


where the contour T° goes from (a — 7-@) to (a + i- ©), a real, a > O in the complex 

s-plane. By a well-known procedure, Eq. (12) leads to 

{Kaaler’) « \ f° o-n{ Kalo’ Ttalo) — Kie)Iafor)\ a 
= , at *an\F)fanlot )\ Oe 


o 


oe ike Gee | [Kio P + rl(oP 








+ Rex Keer)} (12b) 


The integral probably can be evaluated only by numerical methods, and the residues 
of the function (although the singularities, which are in the left-half s-plane, are simple 
poles) can be evaluated only after determining the location and number (which are of 
the order of 2") of poles. The calculation complexity is then clear, even for this particular 
case, r = 1, and an approximate method of solution is sought. 

6. An approximate ‘‘body” solution y'*’ and the complete solution ¢. The comparison 
between the “flat plate’ and “modified flat plate” solutions implies that the main 
contribution to the pressure of the ‘‘body”’ solution will be found near the leading edge- 
body junction, i.e., for small z. With this in mind, an approximation for large s to the 
Green’s function, G.(P; Q), given by Eq. (7e) for the region exterior to the unit circle 
is constructed. Consider a Green’s function of the form 


G(r,6:r’,0’) = K,|s(r? + r” — 2rr’ cos (6 — 6’))'”?] 


7 or 1/2 
+A Kar 4 7a — 77 COs (0 — *’) | (13a) 


where A is to be determined. It is clear, from the invariance properties of the differential 
equation (cf. Sec. 2), that the exact Green’s function cannot be put in this form, since 
the second function on the right side of Eq. (13a) is the fundamental solution (with 
respect to P(r, @)) placed at the inverse point of Q(r’, 6’) with respect to the unit circle. 
But, by the addition formula (ef. Eq. 7b), 


" ] 2r nnd . 8 
K, a(1 + 7a — =7 COS (@ — *’) = K.(sr)I, a 


+2)>> K(s)1.(%) cos [n(@ — 6’)], r> i (13b) 


r 
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Replacing this expression in Eq. (13a) and comparing with the exact Green’s function 
(Eq. 7e), such that the first term is exactly matched, the arbitrary function A becomes 


Pen Is(s) _ Kolsr’) | (130) 


Ki(s) 8 
1( ,) 
7 


Then the approximate Green’s function for the region exterior to the unit circle is 


G.(r, 6:7’, 0’ = K,[s(r" $+ r’? — Qrr’ cos (6 — 6’))' "7 


Ié(s) K,(sr’) ' ] 2r - 
7 Ko} sir 73 cos (6 — 0 (13d 

¥ Ko(s) Io(s/r’) T r r , ) ) 
This approximate Green’s function satisfies the differential equation, Eq. (2a), with 
respect to the point P(r, @), but it is not symmetric in (r, 6; r’, 6’) and does not satisfy 
the boundary conditions. However, comparing term by term, the difference in the 


Green’s functions occurs only under the summation sign, where 


IX(s) 


R K,(sr)K,(sr’) cos [n(@ — 6’)] (13e) 
\n\§) 


in the exact Green’s function (Eq. 7e) has been replaced by 


; / 

Lo(s) I,(s/r’) K,(sr)K,,(sr’) cos [n(@ — 06’)] (13f) 

Ko(s) Lo(s/r ) 

in the approximate Green’s function. Then for fixed 7’ and large s (for Re s > 0) terms 

given by Eq. (13f) approach the exact terms given by Eq. (13e), using the asymptotic 
expansions for the modified Bessel function [11]. 

The approximate Green’s function for the region and boundary for the wing-body 

problem (Fig. 2) can be written immediately 


G(r, 0;r’,0’) = {K,[s(@r? +r” — 2rr’ cos (6 — @’))'”?] 


> 


+ K,[s(r? + r’”? — 2rr’ cos (6 + 6’))'”*]} 


Ii(s) K.(sr‘) J ,- ( 2 l 2r , re, a 
ai a a ed aT rT i 7 COS — 
+ Ki(s) In(s/r’) (Ki) r+ me 2 cos (6 6 )) 


or 1/2 
+ Kr a 7 _ asd cos (6 + a’) } (13g) 
/ 


The approximate Green’s function is a good representation for large s (and r’ > 1); 
this result implies, by a known theorem in Laplace transform theory [13], that the 
solution obtained by using this Green’s function is a good approximation near the leading 
edge-body junction. The transformed approximate “‘body”’ solution is 


os Wa f° 1 I‘(s) K,(sr’) f p ( > l 2r )") 
sy (r,6;3) = — ( + ') ok gee, of Sir” + az — = Cos 
sy” '(r, 85s) } l 78) K's) I9(s/r’) = | ? e 77 COs 6 


é 


: . 1 2r \ ‘as ; 
+ Ko} sir + 7 + 7 COS 8) dr’. (14a) 
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The second term in the braces of Eq. (14a) corresponds to the reflected solution with 
respect to the y, z plane (cf. “flat plate’ solution, Eq. 9c). Now the solution for the 
pressure component, ¢:”’, can be obtained in Region II (ef. Fig. 3) where the approxi- 
mation is best (the reflected solution does not enter here). By the asymptotic expansions 
for the modified Bessel functions 


Ii(s) Ko(sr’) 1 s(r’ — 1)° (1) 
~ Ki(s) Io(s/r’) pa ds {. r it +o s/J° (14b) 








Neglecting the higher order terms 


Ié(s) Ko(sr’) ,, ( _— "] 
K(s) Io(s/r’) Ks r+ 3 — = cos 6 


.* 2 . 1/2 
ww 5 exp - ee )) Viel o(v a z - cos 0) |. (14c) 





r 


For simplicity, the solution is carried out for 6 = 0, i.e., in the plane of the wing only. 
Applying the inverse Laplace transformation to the right side of Eq. (14c) 


£ Mb exp ‘ = a \i.(s ht = : );s} = 0, r’ > (ze —r-+ 2), 
r | ri (144) 


= {lzr’ — (r’ — 1)?)? — (rr’ — 12)7)>7™, l<r’ <(@—r+ 2). 


—_ | 
t 


Then the approximate “‘body”’ solution in Region II is 


so (1 + 3) ler —( — 1° — @r’ — 1)}""* dr’. (15a) 
T rT 


Wa Ji 


As is to be expected, the integral vanishes for (r — z) = 1, z ¥ 0. This is an elliptic- 
type integral which can be expressed in terms of the standard complete and incomplete 
elliptic integrals of the first and third kinds by known methods of reduction [14]. For 
example, in the particular case r = 1, Eq. (15a) reduces to 

gy, (1,0,z) 1 f° -2 21)-1/2 

ee ws | fl + (1 +27) "*}{(1 — (1 + 27)[1 + (1 + 2)7 — 27')} dr. (15b) 


Wa T Jo 
The complete approximate solution, g;'’ + ¢.°’, is shown in Fig. 4 as heavy, solid lines. 
This solution also approaches the correct asymptotic value for large z, i.e., g./Wa = 1. 
No estimate has been made for the error, although in the JPL report such an estimate 
is indicated for a simpler problem, using the same approximate Green’s function for the 
circle (Eq. 13d). Also, additional calculations on the body have not been made; however, 
it is clear that the complete solution approaches the correct asymptotic value very 
rapidly in the downstream direction on the body and in the vicinity of the body on 


the wing.* 
*A comparison with experimental pressure distributions is made in [15] and results, obtained by 


integrating the pressure over finite wings (including cases with swept-back leading edges), have been 


used in [16] on supersonic wing-body lift. 
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An estimate of the complete solution on the body (r = 1) for large z is readily made. 
From Eq. (9b), setting , = 1, the Green’s function for small s becomes 


(K yp! , , \ 
} Ao( sr" ) K,,(sr’) 
< — 2 | 9 en " ~ : eee ~ egies a 3 wy! 9 
= 2 > : cos 2n6? = —s| K,(sr’) 2 
| K,(s) ; K:,,(s) ) + 


\ l 


cos 2n6 


= 16a) 
= Onr'™ ( 


Since K,(sr’) is the dominant term for small s (Re s > 0), the transformed complete 


solution becomes 


sy(1, 63s) 2f : . 
sy = | (1 + 1) (sr’) dr’. (16b) 
r 


Wa wT J; 


Taking the inverse Laplace transform of K,(sr’), the pressure component 


v.(1,0,z) 2 na . : , 9 2 pos 1/2 
7 | (1 + 4 2 — x?) dr’ = 14+ ee ee (16c) 
; 


Wa T Ji T 2 
and 


( 
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AXIALLY SYMMETRICAL JET MIXING OF A COMPRESSIBLE FLUID* 


BY 
&. 1; PAI 
Institute for Fluid Dynamics and Applied Mathematics, University of Maryland 


1. Introduction. The problem of turbulent jet mixing of an incompressible fluid was 
first analyzed successfully by Tollmien' in 1926 by the application of Prandtl’s mixing 
length theory. He solved the following three problems: 


(1) Mixing of a parallel stream with the adjacent fluid at rest. 
(2) Two-dimensional jet from a very narrow opening issuing into medium at rest. 
(3) Axially symmetrical jet escaping from very small opening into medium at rest. 


The corresponding laminar problem has been theoretically analyzed by Schlichting’ 
in 1933. Neither of these solutions, however, holds for finite opening, i.e. at a short 
distance from the opening. 

Kuethe,*® in 1935, extended Tollmien’s results to the case of a two-dimensional jet 
issuing into a medium not at rest and also worked out an approximate method for the 
computation of the velocity profile in the initial part of a round jet issuing into medium 
at rest. Squire and Trouncer,* in 1944, extended Kuethe’s results to the case of a 
round jet issuing into a uniform stream by assuming certain velocity profiles across the 
jet. Abramovich,’ in 1939, extended Tollmien’s problem, i.e. the mixing of a parallel 
stream with the adjacent medium at rest, to the case of a compressible fluid. He 
considered the effects of compressibility due to high temperature and those due to 
high subsonic speed separately and obtained some approximate solutions. Reichardt,° 
in 1941, from the experimental data, suggested the constant exchange coefficient over 
each cross-section of the mixing zone for the free turbulence problem and Gértler,’ in 
1942, reexamined Tollmien’s problems by the application of Reichardt’s theory and 
obtained some improvement in the velocity profiles. 

In 1949, the present author,” investigated the problems of the flow of a two-dimen- 
sional jet from a finite opening of a compressible fluid exhausting into uniform stream 
and of the mixing of two uniform streams of compressible fluid. Both the laminar and 
the turbulent cases were considered. The effects of compressibility due to large tem- 
perature difference and those due to high velocity were considered simultaneously. 

In the present paper, this analysis is extended to the case of an axially symmetrical 
jet of a compressible fluid exhausting into a uniform stream. The flow of the jet is 
assumed to be under full expansion from a nozzle, i.e. the pressure of the flow at the exit 
of the nozzle is exactly equal to that of the surrounding stream. The pressure gradient 
in the jet is assumed to be negligible. Both the laminar and the turbulent cases will 
be considered. 

The first part of this paper is concerned with the laminar flow. The usual assump- 
tions of boundary layer theory adopted to simplify the Navier-Stokes equations. A 
solution by the method of small perturbations is first obtained. -Then the exact solution 
is examined. A numerical integration method is used to compute the velocity and the 
temperature distributions in the jet for the exact solution. 


*Received June 25, 1951. This work was carried out under Contract AF 33(038)-10481 sponsored 
by the Office of Air Research. 
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The second part of the paper is concerned with the turbulent flow. The fundamental 
equation of motion of an axially symmetrical jet is derived by using Taylor’s hypothesis 
concerning the transport of vorticity and Reichardt’s. assumption of free turbulence, 
i.e. the assumption that the exchange coefficient over each cross-section of the mixing 
zone is constant. By suitable transformation of variables, the equations of turbulent 
flow become identical to those of laminar flow. Hence the solution of Part I can be used 
to the turbulent flow provided that the empirical quantity of eddy viscosity has been 


evaluated experimentally. 
THE LAMINAR JET MIXING 


2. Fundamental equations. The equations of motion for the mixing of an axially 
symmetrical jet in a viscous compressible fluid are 





Ou L Ou re) ( an) (1 
= *% PS a r= ) 
P Ox P or Y or 4 or 
Op 
= 0) (2) 
oH] 


These equations are obtained under the usual assumptions of boundary layer theory, 

and the pressure p is assumed to be constant. The x axis is taken along the axis of the 

jet, and r is the radial distance. The quantities u and v are the x and r components of 

the velocity, respectively, p is the density, and y is the coefficient of viscosity of the fluid. 
The equation of continuity is 


O( pu) 1 d(pur) 


: : = 0 (3) 
Ox r or 
The equation of energy is 
fs) = ae 1 a oT du \* 
ou — {C7 ) + ow (CLT) = : (ar ) + (2) (4) 
Ox i or r or or or 


where C, is the specific heat at constant pressure, \ is the coefficient of heat conduction, 
and T' is the temperature. 

If the Prandtl’s number C,u/A is assumed to be unity and C, to be constant, the 
temperature becomes a function of velocity only, i.e. 
u” 
2C 


T=A+Bu- (5) 


» 
where A and B are constants determined by the boundary conditions. The relation (5) 
is well-known for the two-dimensional case,’ but it is interesting to find that it holds 
true also for axially symmetrical flow. 
The pressure is assumed to be constant; hence 
p T, l 
* — = = ~ (6) 
where the subscript 0 represents some reference conditions. 
The coefficient of viscosity may be expressed as 


~J 
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where m may be taken as a constant less than or equal to one depending on the tem- 
perature range investigated. 

3. Solution by the method of small perturbations. In some practical cases, both the 
velocity and the temperature in the jet differ only slightly from those of the surrounding 


stream. It is advisable to find some approximate solution for this special case first. 
Write 
u=UtuUy, v=0,, 
(8) 
P= p tp, B= bo +m 
where u bsg Oy Ry yy Ke Oe es KE ee 
Substituting Eq. (8) into Eq. (1) neglecting the higher order terms, one has 
ou »loa Ou 
- ‘=a’ — (; ma) (9) 
OX ror or 
W here 
= Ho 
Uo Po 
The initial conditions at the exit of the nozzle may be assumed as follows: 
z= 0 
Uu, = Uy = constant, o<r <= ¥F, (10) 
u, = 0, oe a 
The solution of Eq. (9) with the boundary conditions of Hq. (10) is” 
Uy ma : , 
4 =n] e Jo(dr)J (Aro) ad, (11) 
Uro Jo 


where J, and J, are the Bessel functions of zero and first order respectively. The results 
are plotted in Fig. 1. 
Atz = 0 
oa lL. fp ot, 
= 1 | Jo(dr)Ji(dra) dd = (12) 
~ ie 0, <r. 


This is the well-known discontinuous integral of Weber and Schafheitlin'’. Thus the 
boundary condition (10) is satisfied. 

4. Exact solution. To solve the problem more rigorously, one has to resort to Eqs. (1) 
and (3). By introducing the stream function y, which is defined by 


Oy _ oy 


,e err, Bae ee (13) 


or 

where p* = p/p, , u* = u/u , 7* = r/ro and z* = 2x/ro ; the equation of continuity, 
Eq. (3), is automatically satisfied. 

If the independent variables x and r are changed to z{ and y. [Equation (1) becomes 


au* 0 2 ou* 
aa = ay (ur tutes? au") (14) 
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Fig. 1. Velocity distributions in axially symmetrical jet by small perturbation solution. 
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where x* = a’ x*. The initial conditions are now 
a =e. T* = THY), u* = us(y) (15) 


Equation (14) can be integrated numerically by a finite difference method starting 
from the initial conditions (15) with the help of (5), (6) and (7). In the numerical 
integration, besides (15), the condition on the axis of the jet must be known which is 
found to be: 

ou* Qu* du* 


2 (16) 


ox* p*u* or 
As soon as u* and p* are found in terms of y at given x*, the radial distance r* can be 
computed by the following formula 
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A typical example has been set up for IBM Machine calculation by Dr. H. Polachek 
and Mr. T. 8. Walton of the Naval Ordnance Laboratory. The boundary condition at 


the exit of nozzle (2* = 0) are 


. u* = 1.0 


I” = ae < is TT? = 10 79> 1 


* 


u 


Ba 


p* = 0.5 p* = 10 
The following relations were also used in the numerical example: 
T* = —9.66 + 11.26u* — 0.60u*’, p* = 7" 


This represents a hot jet of flow in a cold stream. The result is shown in Fig. 2. To 
compare with the small perturbation solution, it shows that the velocity on the axis of 
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Fic. 2. Velocity distributions in axially symmetrical jet by numerical integration. 
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the jet decreases more rapidly for the case of the hot jet. These results is qualitatively 
the same as that in case of two-dimensional jet flow.* 
THE TURBULENT JET MIXING 


5. Fundamental equations for turbulent jet mixing. The differential equation of mo- 
tion for steady axially symmetrical flow of free turbulence, where viscous force can be 
completely neglected, is 


Ou Ov Op 
us = (18) 
“ Ox re or Ox 
If we put 
U=Un tu, v= v0, + v’, 
(19) 


P= pm t p’, P=PDrntp’, 


where u, v, p, p are the instantaneous values of the x and r components of velocity, 
pressure, and density, respectively, u,, , Un, Pm, Pm , the corresponding mean values and 
u’, v’, p’, p’ the corresponding fluctuating values. 

Substituting (19) into (18) and averaging it, one has 


‘au’ ou. Ou 
Palln| ‘) + ake | + | (p u')a( 24) 
Ox = \ or) = \Or = 
, Ou’ , ou’ . Ou 
+ ual p ) + pa(u : | + co (2 
OL I « OX / m OT? as 
4 \ wd 4 r\ 
, Ou , OU u 
+ o{ p - ) + pal v : ) |+ | (ow = ) 
. ov. OT sa OL / » 
, ou’ ‘Op 
— (o’ au | | a P) ; (20) 
\ m lm 


o7 \Ox 


where (__),, denotes the mean value of the quantity in the parenthesis. 
By comparison of the order of magnitude and using the ordinary boundary layer 
assumptions one has 


au Ou er. , ou’ ap 
pues ) +" pat a( 2) + (pv (2) + pa(v = ) _ -(% (21) 
\OX/ m a ee Or] m nT la OX 


im 


In order to find the relation between du’/dr and the mean flow, we use Taylor’s 
modified vorticity transport theory,’’ i.e. 


Du’ lod ‘a 
— = Aw = 1-- (2) (22) 
or r or Or] m 


where Aw is the change of vorticity of mean flow in r direction and I is the mixing length. 


Furthermore we put 
dp 
= (2) : 
p or], (23) 


tw 
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Substituting (22) and (23) into (21) and put (dp/dz),, = 0, one has 


au -4 eres Pes 
pata) T Pmt { me) ~ <> ar Ee | (24) 


where « = —(lv’),, is the exchange coefficient of the turbulent flow. According to Reich- 
ardt’s assumption, we assume that « is independent of r in free turbulence problem. Hence 


ou ou 1 oa Ou a 
patra) TF Pmt { x) 7 or | eonr(24) | (25) 


From (25) we see that the turbulent shearing stress of compressible fluid in the present 





problem may be defined as 


Ju . 
Tor = ev.(2 ) = —p,,(U’'v’) m (26) 


_ 2) 
un’ = (2 - 


The expression of (26) is in agreement to that obtained by Frankl,’* but we use 
different assumptions. In Frankl’s analysis, he assumed that p’ is small by comparison 
with p,, which is not a good assumption in high speed flow. 

The energy equation for free turbulence case is (from hereon we drop the subscript 


where 


m for mean values) 


a(C,T) ac,T) 12a ( a) (2) - 

uy 2 py = —-— IC 20 | €p\ — (27 

4 Ox T Pp or r or oP ar oe or ) 

According to mixing length theory, the Prandtl] number for free turbulence is equal to 
unity, equation (5) holds true for turbulent flow too. 

6. Solutions of the equations. According to similar arguments for two-dimensional 


8 . 
case, we may write 

z n 

ro 


where ¢€) is a constant to be determined by experiments. 
Equation (25) becomes 


Ou ou 1 od ( 2) 
— y— = e z*" — — ,— (29) 
= Ox pt or - r or a or 


Under the assumption of small perturbation equation (29) becomes 


dus — gor 2 (, 20) (30) 


ax ti‘ ror\ or 


0 


In terms of variables z* and y, equation (29) becomes 


ou* _ €0 *” 


= 2 ( a oe 2u*) (31) 
az* ~ uty ay \? a 3 
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Now we introduce the new independent variable 


(n+1) 
"ig 
n+ 1 
for x*, equations (30) and (31) become respectively 


Ou , bw Ou : 
es € z (r= ) (33) 


X 


OX ~~ Ulo r* Or* or* 
and 
= , i 
ou a 0 aoa et CU (34 
te — ieee 34) 
OX Uol'yo OW or / 


are identical to equations (9) and (14) respectively, provided 


Equations (33) and (34) 
Therefore the solutions of the laminar 


that ¢, is used for uo/p) , X for 2* and p* for p*. 
jet mixing can be applied to the problem of the turbulent jet mixing of the same boundary 


conditions, provided that proper characteristic constants €) and n are chosen. 
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TORSION OF A CIRCULAR CYLINDER HAVING A SPHERICAL CAVITY* 


BY 
CHIH-BING LING 


Aeronautical Research Laboratory, Taiwan, China 


Abstract. This paper presents a torsion solution, based on the Michell-Féppl theory, 
for an infinite circular cylinder having a symmetrically-located spherical cavity. Nu- 
merical values are also given to show, in particular, the effect of the cavity on the 
maximum shear stress in the cylinder. 

The general theory. Let (r, 6, z) be the cylindrical coordinates of a point. For con- 
venience, 7 and z will be considered as dimensionless quantities each measured by a 
unit of certain typical length a. 

The solution of the torsion problem for a solid of revolution, requires a function y 
which satisfies the following differential equation.’ 

of. 2 4 ee wt (1) 
or” r or Oz 
in which the axis of revolution is taken as the z axis. The two non-vanishing stress 


components are expressed by 


J d 
At a (2) 


a Ta —— ae _ 
where uw is the modulus of rigidity. The condition that the surface is free from traction 
takes the form 
wv = const. (3) 
on the entire surface. 
The constant for ¥ at the internal surface will henceforth be taken as zero as no 
generality is lost in doing so. The terminal couple acting on the cylinder is then equal to 


T = 2Qnrua*y, , }) 


where y, is the value of ¥ at the external surface. 
The above theory is due to Michell and was rediscovered by Féppl. 
Solutions in cylindrical and spherical coordinates. If we put 
y= ry (5) 


where WV is a function of r and z, the preceding differential equation becomes 





yy 1 ON 4y ry 
e 4 ad «a A Co a (6) 
or” r or i Oz 
Let us assume a normal solution of V in the form 
sin kz 
vY = F(r) ‘ (7) 
cos kz 


*Received April 11, 1951. 
iSee A. E. H. Love, Mathematical theory of elasticity, 4th edition. Dover Publications, New York, 


1944, pp. 325-326. 
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It then appears that the function F satisfies Bessel’s equation 
dF , 1dF at 
s+ —ik + 3/F = 0. (8) 
dr r dr r 
This equation is satisfied by J.(kr) and A.(kr) which are modified Bessel functions of 
order two, of the first and second kinds, respectively. 
The above solution holds for any value of k. A more general solution of y is therefore 
expressed by the following four integrals: 
“e I.(kr) sin kz 
yar | f(k) x dk, (9) 
es K.(kr) cos kz 
where in each integral f is an arbitrary function of k. 
Now we proceed to find the solution in the spherical coordinates (p, ¢, 9) which are 
connected with the cylindrical coordinates by 
z= pcos®d, r= psin®¢, 6 = 80. (10) 


Consequently, the differential operators are transformed to 


6] 0 sing 0 
-_ =~ COSB@Q~ oa a. 
Oz 0p p odo 
(11) 
0 ; 0 , cos¢d ~A 
~— =snqQo- — ee 
or 0p p Op 
The differential equation for Y becomes 
ov 2 ow 1 ow cot d dV hy 
—_ F “ss iT oo = AE ) (12) 
Op p Op p oo p fefe0) p sin @ 
Assume a normal solution in the form 
VY = G(¢)p’. (13) 
It follows that the function G satisfies the following equation 
&G dG . oc 
s+ cotd + {n(n + 1) — =3— 76 ='0. (14) 
dd do q sin’ @) 


This equation is known as Legendre’s associated equation and is satisfied by P,, (cos ¢) 
and Q? (cos ¢), which are Legendre’s associated functions of degree n and order two, of 
the first and second kinds, respectively. 

It is noted that the same Legendre’s associated equation is obtained if we assume 
instead 

G(¢) 
= —, 

p 


Since the preceding solution is valid for any integral value of n including zero, a 


-_ 


V (15) 


more general solution of y is expressed by the following four series: 


" P:(cos ¢) p- 
y=rDA x ; (16) 
sie Q;(cos d) p 


where A, are parametric coefficients. Note that P; vanishes when n = 0 or 1. 
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Method of solution. Now, consider an infinite circular cylinder of radius a (i.e., 
r = 1) having a symmetrically-located spherical cavity of radius Xa (i.e., p = A), which 
is under torsion by terminal couples about the axis of the cylinder. For convenience, 
the centre of the cavity will be placed at the origin. 
The method of solution is to assume ¥ as being composed of two parts as follows: 


YVrwrh, (17) 


where y, is the solution of a corresponding cylinder without a cavity, while y, is an 
auxiliary solution which vanishes on the surface of the cylinder. The latter is added to 
yY) so that the remaining boundary condition on the surface of the cavity is satisfied 


by adjusting the parametric coefficients involved in the function. The function yp is 
vy = Arar’, (18) 


where 7 is a constant. Since ¥, vanishes on the surface of the cylinder, it follows that 
the terminal couple is equal to }aura*. 

The auxiliary function y, may be constructed by combining linearly a set of functions 
each of which satisfies the given differential equation and vanishes on the surface of 
the cylinder. It appears that the set of functions can be derived from the following 
particular solution for y, namely: 


y* = 2 P3(eos ¢) + r° | f(k)I (kr) cos kz dk. (19) 

It is obvious that the function so constructed satisfies the given differential equation. 

Note that this function is even in z or cos @ and, besides, it has a singularity at the origin. 

According to the definition by Hobson,’ ?3 (cos ¢) is equal to 3 sin® ¢. Hence ¥* 
vanishes when r = 1 provided that 


| f(k)I.(k) cos kz dk = — ? ee (20) 
Jo (1 + 2)"" 
A Fourier transform gives 
2: “= cos kz dz “Kk 
oe 2 3 _cos k td = _ 2k’ K,(k) (21) 
wIA(k)Jo (1 +27)” wl ,(k) 


The last result is a particular case of the integral considered by Poisson and 
Malmstén.* The function ¥* is thus fully determined. It is obvious that differentiation 
with respect to z gives functions with the desired properties on the surface of the cylinder, 
but odd derivatives must be excluded since they are not even in z and cannot enter 
into the required solution. The function y* itself may also be included. The set of func- 


tions is therefore 


2 pe 4s 
y*, hs 0 : peek, (22) 
Oz Oz 
*See E. W. Hobson, Theory of spherical and ellipsoidal harmonics, Cambridge University Press, 1931, 
p. 94. Also ef. E. T. Whittaker and G. N. Watson, Modern analysis, Cambridge University Press, 1927, 
p. 325. The definition given in the latter is slightly different though it is also attributed to Hobson. 
See G. N. Watson, Theory of Bessel functions, 2nd edition, Cambridge University Press, 1944, p. 185. 








152 CHIH-BING LING [Vol. X, No. 2 


Since the use of any constant multiplier does not affect the desired properties, we 
may write the set of functions as follows: 
l a" y* 
73, = : v (23) 


(2s — 2)! a, 


the initial function being y% = y* 
To express the functions in terms of the spherical coordinates, the following relations 


> 4 
are useful. 


d ry =r sa (2s —2)!r?_, 
; Px(cos o) (7 = 3r vert tert in P;,(cos ¢) (24) 
Oz pP ) Oz ~ \p p 
and 
l i 1 (kp)” 22 95 
bAKT) C08 he = = 22 (4) ; P2,( COs ¢). (25) 
: (Qn + 2)! 


The latter expansion is valid in the neighborhood of the origin. Hence, it follows that 


for s > l, 


P2 (cos $ ; " ; = 
yar ~ 7? FP? (cos ¢)p", 26) 
p- 1 
where 
rs {2n + ~ o o7 
Yo, = (- 2 
\2s—2/9 Pe 
in which, for s > 4, 
22 ; K.(k) 
, = 2 dk. 28) 
x s!. [,(k 


Note that the asymptotic value of o, ise jual to unity. 





By means of the above set of functions we construct 


¥, = a 2, AsV. , (29) 


where A,, are parametric coefficients to be determined; the factor ra being introduced 


to render A,, dimensionless. The function y is then equal to 


rar Lf A y 2 . 
y = 19 P>(cos ¢)p + rar 2 a be Qo,Ao \P2.(cos od). 30) 
2 p 
Hence the boundary condition that y = 0 when p = X is satisfied if, for n > 1, 
| . A on WY 2 
N’6 =~" > “a,,Az, = 0 
12 d 
or (31) 
| int] 2 
| = = 12 A 0 — nN - a he ; 
where 6, = 1 or 0, according as m = n, or m # Nn. 


‘For the first relation see E. W. Hobson, Loe. cit., p. 105, formula 36. For the second relation, ef. 
T. M. MacRobert, Spherical harmonics, 2nd edition, Methuen, 1947, p. 362, Ex. 63. Also see the recur- 


rence formulas for Legendre and Bessel functions. 
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The above system of linear equations may be solved by successive approximations 
as follows. Write, for n > 1, 


Am (32) 


where 


r) 


AS = TM bia, AS? =O DE May A (33) 

The validity of the above solution naturally depends upon the convergence of the 

series (32). To establish convergence, an inequality for the coefficients is needed. How- 

ever, for the sake of brevity, no proof will be given here. From physical considerations 

alone, it seems likely that there will be convergence as long as the boundary of the 
cavity does not touch the surface of the cylinder or \ < 1. 

Effect of cavity on angle of twist. To investigate the effect of the cavity on the angle 


ol twist, we have for ¥ 


1 dw, 1 dy 
a z- om 4 
a Oz ar or : (34) 


where w, is the angle of twist due to y alone. This implies that 7 represents the angle 


of twist per unit length of a circular cylinder without a cavity. While for y, , we have 


1 oy _ x - 5D Ay, He, (35) 


a Oz ar” = J 


where w, is the angle of twist due to y, . Integration gives 


ai >> Az S / vi, dz. (36) 


It appears that there is no contribution from y* except the integral part in ¥% . The 


latter gives 


iraAj od . ff f[- k* K(k) I2(kr) oe 

wo = ta | - Lb cos kz dz dk 
47a Az _. * k* K(k) (kr) sin kz > 
ee oe I “I(k) -” 


(37) 


_ _2raA, Li Sk? K(k) 1 (kr) | 
= A Lim Taf 


= —1l67raA,. 
i.e., owing to the presence of the cavity the angle of twist of the cylinder is increased 
by an amount —167aA, : 

Numerical examples. Numerical examples will be given for several values of \. Values 
of o, are given in Table I, which are computed by numerical integration with the aid 
of Gregory’s formula.® The coefficients ““a., are then readily obtained and shown in 
Table i. 


. E. Whittaker and G. Robinson, Calculus of observations, 4th edition, Blackie 1948, p. 143. 


we 
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Values of A,, found by successive approximations are shown in Table III. The stress 
at any point in the cylinder can now be readily computed. In particular, the shear 


stress across the minimum section z = 0 is given by 


[Qn — 1 (Qn + 2r""|_. en 
To, = —uta >, A,,§———— + ae P?,(0), (38) 
n=l c.rlC } ° j 


where 


(2n + 2)! 


* (0) = (-1)"' = : 
: 2°"(n — 1)!(n + 1)! 


(39) 


This stress is shown in Table IV and also graphically in Figure 1, where the known 
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Fic. 1. Shear stress across minimum section 


result for \ = 0 is also included. It appears that in the cylinder under consideration 
the maximum shear stress always occurs on the external boundary across the minimum 
section. Now, if we define the stress concentration factor AK as the ratio of the maximum 
shear stress to the constant shear stress across the same section which would give the 
same couple about the axis of the cylinder, i.e., 

4(1 — 2°) 


A = —— (Fea (40) 
dura 


then the following results are obtained. 


2a a a ee 
K | 1.383, 1.3827, 1.292, 1.220, 1.139, 1 
The value of K in the limiting case \ = 1 can readily be visualized from physical 


consideration of the cylinder. Figure 2 shows the results graphically. 
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Fic. 2. The stress concentration factor K versus X. 
TaB_e I. CoEerriciENT o; 
} 6 s 10 12 
Os 12.81758 $.15488 2.75876 2.20244 1.91108 
14 16 18 20 22 24 
Os 1.73288 1.61286 1.52653 1.46174 1.41078 1.37056 
TaBie II. CoerriciENT *"ao, 
2. 2n = 2 2n = 4 2n = 6 2n = 8 2n = 10 2n = 12 
2 8.011010 -6.4920 « 107? 1.0776 107? | —2.1508 10-3 4.6657 X10-* | —1.0577 X10~* 
4 9.7380 «107! 2.9080 X 107! | —9.6787 X10? 3.079410? | —9.6248 X 10-3 | 2.9532 x 10-3 
6 7.5485 «107 | —4.5167X107 2.3095 X 107! | —1.0587 K10™ 4.4791 X1072 | —1.7819 X10? 
8 +.5167 107 1.3111X10" | —3.1762X<10" 1.9708 X10™! | —1.0810X10~ | 5.4032 x 107? 
10 2.3095 x 107! -3.1762X<10™ 3.1673 X10" | —2.5481 K107 | 1.7561 X10— | —1.0756 X107 
12 1.0587 X 107! 1.9708 X10" | —2.5481«K10™ 2.5755 107! | —2.1750 107 | 1.6022 107! 
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0718 
9045 
8434 
3619 
2796 
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TaBLeE III. Coerricient 

X 1 
10~° 3.4407 X 10 2 
10714 1.13848 * 10 3 
10 2.3256 107! 3 
10 5.7307 X 107~% { 
10~-? 1.53847 X 107" 5 
10-9 $.2955 * 10°29 Ss 
TABLE [V. SHEAR STRESS 79;/p7a 
0 r I r 1 
000 0.208 0.418 
250 0.377 0.518 
500 0.584 0.674 
790 0.792 0.858 
000 1.000 1.006 
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5176 
3912 
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Ac ROSS 


10 

10-7 
| ld 
10-™ 
10-3 
10-* 


0.643 
0.696 
0.799 
0.919 
1.046 
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2.0889 * 10 
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SOIL MECHANICS AND PLASTIC ANALYSIS OR LIMIT DESIGN* 


BY 
D. C. DRUCKER anp W. PRAGER 


Brown University 


1. Introduction. Problems of soil mechanics involving stability of slopes, bearing 
capacity of foundation slabs and pressures on retaining walls are often treated as problems 
of plasticity. The soil is replaced by an idealized material which behaves elastically up 
to some state of stress at which slip or yielding occurs. The shear stress required for 
simple slip is often considered to depend upon the cohesion and linearly upon the normal 
pressure on the slip surface. In more complete plane investigations an extended Coulomb’s 


rule is used,** (Fig. 1) 


Co ys 
R=ccos¢g — > sin ¢, (1) 
where F is the radius of Mohr’s circle at slip, the maximum shearing stress [(¢, — a,)°/4 + 
r.,|°; ¢ is the cohesion; ¢ cos ¢ is the radius of Mohr’s circle at slip when the mean 


normal stress, (¢, + o,)/2, in the plane is zero, and ¢ is the angle between the tangent 
to the Mohr’s circles at slip and the negative o axis. 


§ 














Ty toy 








Fic. 1. Mohr-Coulomb hypothesis. 


The purpose of this note is to discuss the implications of assuming the soil to be a 
perfectly plastic body. The validity of this assumption is not at issue. No account is 
taken of such important practical matters as the effect of water in the soil or of the 
essentially different behavior of various constituents such as clay and sand. All that is 
sought is a theory consistent with the basic assumption. Comparison of the predictions 
of such a theory with the actual behavior of soil will give an indication of the value of 
the idealization. 

*Received Nov. 19, 1951. The results presented in this paper were obtained in the course of research 


sponsored by the Office of Naval Research under Contract N7onr-35801 with Brown University. 
**KX. Terzaghi, Theoretical soil mechanics, John Wiley and Sons, 1943, pp. 22, 59. 
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2. Yield function and stress-strain relation. A yield function which is a proper 
generalization of the Mohr-Coulomb hypothesis (1) is 
f = aJ, + J," = k, (2) 
where a and k are positive constants at each point of the material; J, is the sum of the 
principal stresses: 


J, =o, +o.+o, =0o,+06,+ 0 O11 + Foo + 633 C, 


oa = 58 lio, — o,) + (¢, —¢.) + (6. —o,) | +r,+¢7.4+7 
The stress deviation is defined as s;; = o;; — (J,/3) 6;; , where 6,; is the Kronecker delta, 
zero fori # j and unity for? = j. For example, s,, s, = 3[o, — 3(o, + ¢.)], Si2 
Gae = Fay 
As is well known, the yield surface, f = k, in principal stress space is a right circular 
cylinder for the Mises criterion (a2 = 0). For a > 0, the surface is a right circular cone 


with its axis equally inclined to the coordinate axes and with its apex in the tension 
octant. A generalization of (1) as a modified Tresca or maximum shear stress criterion 
is also permissible, but will not be discussed here. It leads to a pyramid instead of a cone. 

According to the concept of plastic potential,* the stress-strain relation corresponding 


to the yield function (2) is 


ej; = A Of/de;; , (3) 
where ¢/” is the plastic strain rate and X is a positive factor of proportionality which 
may assume different values for different particles. Substituting the expression (2) for 


f, we obtain 


ef} = \lad,;; + 8;;/2J.2'°]. (4) 
A very important feature of Eq. (4) is that the plastic rate of cubical dilatation is 
é;; = Baa. (5) 


Equation (5) shows that plastic deformation must be accompanied by an increase in volume 
if a ~ 0. This property is known as dilatancy.** 

3. Collapse or limit design theorems. It has been shown previously*** that when 
the limit load is reached for any body or assemblage of bodies of perfectly plastic ma- 
terial, collapse takes place at constant stress. This means that at the instant of collapse 
the strain rates are purely plastic. The assumption is made, just as in most problems 
of elasticity, that changes in geometry are negligible. 

As we shall deal only with collapse or limiting states, e/? equals e!; , the total strain 
rate. 

The major limit theorems for stress boundary value problems are: 

(7) collapse will not occur if any state of stress can be found which satisfies the 


*R. v. Mises, Mechanik der plastischen Formaenderung von Kristallen, Z. angew. Math. Mech. 8, 


. 
161-185 (1928). 
**\I. Reiner, A mathematical theory of dilatancy, Am. Journ. of Math. 67, 350-362 (1945 
***T) C. Drucker, H. J. Greenberg, W. Prager, Extended limit design theorems for continuous media 


Q. Appl. Math. 9, 381-389 (1952). 
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equations of equilibrium and the boundary conditions on stress and for which f < k 
everywhere; and 

(77) collapse must occur if for any compatible flow pattern, considered as plastic 
only, the rate at which the external forces do work on the body equals or exceeds the 
rate of internal dissipation. 

The theorems are valid when frictional surface tractions are present if there is no 
slip or if the frictional forces are known. 

Discontinuous stress or flow states are permissible and are generally convenient for 
computational purposes. 

4. Plane strain. The yield function (2) will first be shown to reduce to the Mohr- 
Coulomb rule in the case of plane strain. In this case €{; , €{3 , €53; (Or €2 , ¥22 5 Yj2 iN en- 
gineering notation) vanish. Therefore, from (4), s;; and 823(7,, , T,:) are zero, and 


833 = —2aJd,’’. (6) 
Thus, 
J, = (O14 + O22) — 3aJ 5 (7) 


J, = | («: = 2) + 2,|/ a — 3a’). (8) 


Substituting (7) and (8) into the yield function (2), we find 


and 


f =3a% = oo 4 (1 — Sa) J? = k 
I 3a , + «, | (22 _ es) , | . 
;= ——n173 <_< 5 9 
(i — 3a) -~G@—3e)" 2 + 2 T tr ) 


Equation (9) becomes identical with Eq. (1) if we set 


k 3a ‘ 
(Il — 12a”)! ia (i- 3a) ae (10) 


and hence 


(1 — 12a’)'” 
q— 3a*)'7 = COS ¢. 


A clearer picture of the meaning of the plasticity relations is obtained by considering 
a plane velocity pattern as shown in Fig. 2. The upper portion is supposed to translate 


t 
a 





Fic. 2. “Simple” slip. 
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as a rigid body while the lower part is fixed. If the transition layer between is made 


thin enough, Fig. 2 represents simple slip. The essential feature is that «f = 0. As €! is 
zero, €, cannot be zero. In fact, from Eq. (5), ef = 3a and the upper shaded region 


not only moves to the right but it also must move upward. Such volume expansion ac- 
companying slip is a well known property of reasonably packed granular material. It 
is an interesting and in fact a necessary result of plasticity theory that if the yield point 
in shear depends upon the mean normal stress, slip is accompanied by volume change. 
This necessary relation is an essential difference between a plastic material and an 
assemblage of two solid bodies with a sliding friction contact.* 

Figure 1 shows clearly that the shearing stress on the slip surface is not the maximum 
shearing stress R but is R cos ¢. 

It is interesting to observe that a is not completely free because 12a° cannot exceed 
unity. As shown by Eqs. (10) this is equivalent to restricting sin ¢ to be no more than 
unity. 

5. Rate of dissipation of energy. The rate of dissipation of mechanical energy per 
unit volume is 

Obes SO 55K a) = Xt = XE. (11) 
Oo; ‘ 
As the dissipation is to be computed from a plastic strain rate pattern, it is desirable 


to express \ in terms of strain rate. It follows from (3) that 


hs , Of af . ae 1 
€; =X = ——. (3a +; (12) 
O0;; Oo; 2 
Substituting in (11) 
k(e,,</,)'"" c COs ¢ ” ‘ , : ~ pA 
D = = ——5 7 [Qe + 2 4+ 22° + ye + yt (13) 
(3a° + 1/2)” (1 + sin® yg)?" ‘ : me Ye] 


which for the plane strain case may be written as 


sii (14) 





D =c cos gl (e — e)° + vy 
For purposes of calculation it is convenient to consider a transition layer as in Fig. 2 
to be a simple discontinuity. However, for a ¥ 0a discontinuity 6’ in tangential velocity 
must be accompanied by a separation or discontinuity 6v’ in normal velocity. A transition 
layer of appreciable thickness therefore must be present in a soil while there is no need 
for such a layer in a Prandtl-Reuss material where a = 0. 
The rate of dissipation of work per unit area of discontinuity surface is 


D, = réu’ + civ’ (15) 


where 7 is the shearing stress and o the normal stress (tension is positive) on the surface. 
The same result can be obtained by taking the limit of tD, Eqs. (13, 14) as ¢, the thickness 
of the transition layer goes to zero while fe{/, = thy!, = }6u’ and tej, = te) = dv’. It 
follows that 

, ‘ 

- §u’ = = sa73 Ou’ = Ju’ tan ¢. (16) 
1 (1 — 122°) 





*D. C. Drucker, Some implications of work hardening and ideal plasticity, Q. Appl. Math. 7, 411-418 
(1950). 
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Also 
ku’ 
DD, = mis = ¢ Su’. (17 
(1 — 12a°)'”” 
These results apply equally well to a curved surface of discontinuity because there 
is no essential distinction locally as the thickness of the transition layer approaches zero. 
6. Dilatation. The rate of cubical dilatation, e/; , which accompanies plastic deforma- 
tion is given by Eq. (5). Substitution of (12) leads to 


3a 


, sr \i/2 
Ge = DT nis (Qe i€ 
(1 + 6a°)'” iets) 
(18) 
- : “ ‘4 172 (Qe? ,¢/;)' : 
(1 + sin’ ¢) 
lor plane strain, we obtain the rate of dilatation 
, , sin ¢ 2 2 2\1/2 
V=e+e= a 175 (Des + el” + y2s)' (19) 
(1 + sin’ ¢) 
or 
A’ = sin gl(e. — &)? + yij'? = I’ sing, (20) 
where I’ is the maximum shear rate. 
In terms of principal strain rates for plane strain, e{ , €: , Where ef > € 
: 1+sing , ¢ 
ef = ef ——— = ei tan’ (45 + 5]. (21) 
l1—sing 2 


7. Simple discontinuous solutions. Any soil mass, in particular the slope shown in 
Fig. 3, may fail by rigid body “sliding” motion in two simple ways. As is well known, 
the surfaces of discontinuity, idealizations of discontinuity layers, are planes and circular 


cylinders for a Prandtl-Reuss material, a = 0. When a # 0, a plane discontinuity 


surface is still permissible but the circle is replaced by a logarithmic spiral which is at an 


CENTER OF ROTATION 

















Fic. 3. Rigid body “slide”? motions. 


angle ¢ to the radius from the center of rotation. The circle is not a permissible surface 
for rigid body “‘sliding’”’ because a discontinuity in tangential velocity requires a separa- 
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tion velocity. For a “‘slide’’ discontinuity the angle between the velocity vector and the 
discontinuity surface is ¢, as is seen from Eq. (16). 

8. Example: Critical height of vertical bank. The computation of the critical height 
H.. of a vertical bank, Fig. 4, will serve as an illustration of the procedure consistent 
with plasticity theory. No attempt will be made here to get an exact plasticity solution 
to this problem for even if it could be obtained it would not be helpful for more compli- 
cated slopes. Instead, the upper and lower bound techniques of limit design will be used. 











x H tan fh 
bu’ 
+ H ” 
oy = -Wwy y 

Sy = Ty =O bv’ 

— SS ae we a 
6y=—-wy 
Y dv’ bu’ 
Cy |6x i w(H-y) 

anna” 9 
Fig. 4. (a) An equilibrium solution, (b) Plane, (c) Log. spiral compatible solutions 


Denoting the specific weight of the soil by w, Fig. 4a shows a discontinuous equilib- 
rium solution in which the maximum shearing stress is wH/2 at the lower ground level, 
where the average normal stress is —wH/2. It should be remembered that the stress 
field need bear no resemblance to the actual state of stress. Equilibrium alone must 
be satisfied in addition to the yield condition (1) or (2) for H to be safe. From Eq. (1) 


wH _ wH 


as C COS & + 9 sin g 


or a lower bound on the critical height H, is given by 
2c COS @ 2c re) 
“ S¢ “ ~ BY; 
H,.>H = -— = —tan (45+ = (22) 
wil—sing Ww 2 
An upper bound on H, may be found by taking a plane slide as a velocity pattern, 
Fig. 4b. Equating external rate of work to the dissipation Eq. (18) gives 


wH*tanB 6u’ ae 
, — cos (¢ + B) = bu’, 
2 COS ~ cos 8 
so that 
os COS ~ 
|; i (23) 


w sin B cos (yg + 8) 
minimizing the right hand side gives 


B=45-—£ (24) 
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and 
te ¢) 
a ‘ 5 25 
nm. <= tan (454 2 (25) 


The usual theory furnishes the right-hand side of (25) as the critical height. However 
there is a factor of 2 between the upper bound (25) and the lower bound (22) for H, . 

The upper bound can be improved by considering a rotational discontinuity (log- 
arithmic spiral, Fig. 4c) instead of the translational type. 

9. A modified stress criterion. The large discrepancy between the lower and upper 
bounds derived from reasonable stress and velocity fields (Fig. 4a and Figs. 4b, c) has a 
simple explanation. An improvement in the static or equilibrium answer would require 
adding stresses which will result in tensile stress in the soil. Nothing in the Mohr- 
Coulomb criterion, Fig. 1 rules out a moderate amount of tension. If it is felt that soil 
cannot take any tension, then the yield criterion must be altered to take into account 


Co + oO > (2 7 *,) ” 2] : on R (26) 


\s before, if rigid body discontinuous solutions are considered, one strain rate, say 


explicitly that 


bo 


, 


e{, , is zero. The relation between the normal and the shearing stress on a surface of 
discontinuity is now given by Fig. 5 instead of Fig. 1. The yield function and the re- 
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| ¢ > 
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Mh a al 

ges y 
Fic. 5. Modified criterion for stress on discontinuity surface Ro = c tan (45 + ¢/2) 


sulting velocity relations previously derived apply only for the straight line portion of 
the limiting curve. 

In general, the rate of dissipation, D, , given by Eq. (15), may be interpreted as 
the dot product of a stress vector having rectangular Cartesian components o, 7 with 
a velocity vector having components 6v’, du’. The stress vector is shown on Fig. 6 as 
OP or OQ and the velocity vector is normal to the limiting curve. The normality condition 
holds because o and 7 are the only stress components which do work.* Considering the 
stress vector OQ to be the sum of OO’ and O’Q, we have 


f , tan(45 + ¢/2) ae 
ee (27) 


_= >.< | e,,/\2 / J ad | iene ah — = 
D I | LCu’) + (év’)"] 6 { céu tan(45 ++ ' 


> 


- . * ° . &F 
*D. C. Drucker, A more fundamental approach to plastic stress-strain relations, to appear in_Proceed- 
ings, Ist U. S. National Congress for Applied Mechanics. 
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where 


If Fig. 


continuity line as 
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dv’ 
du’ 


Pp 








6. Rate of dissipation on discontinuity surface is adv’ + réu 


év’ = 9e 
tan 6 = — 2? tang (23) 
ou 


c tan (45 + 9/2) (29) 


R, 


tb is modified by taking the angle between the velocity vector and the dis- 


@ instead of ¢g, the results are not altered. However, if a very crude 


Jicture is employed of rotation about a point below the lower ground level, Fig. 7, the 
; I} £ 


upper bound is ea 


sily reduced to 








I 5 ; - 
o. =< = : c tan (45 -+- ¢) (30) 
7 
J) 
—Vy-- O~0 
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a | 
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Rotation about P, soil of zero tensile strength. 


Fic. 7. 
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Further reduction is possible with other discontinuous patterns but the concept of 
a soil unable to take tension clearly requires much additional study. 

10. Conclusion. It has been shown that the implications of plasticity theory and 
limit design for soil mechanics are indeed far reaching. Circular sliding surfaces are re- 
placed by logarithmic spirals in a material whose yield stress in shear depends upon 
normal stress. Volume expansion is seen to be a necessary accompaniment to shearing 
deformation. Low critical heights are found for slopes when the soil is assumed to be 


unable to take tension. 
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A PROBLEM OF FINITE BENDING OF CIRCULAR RING PLATES* 


BY 
ERIC REISSNER 


Massachusetts Institute of Technology 


1. Introduction. We consider axi-symmetrical transverse bending of a circular ring 
plate under the action of transverse edge forces. Let a be the inner radius of the plate 
and a + b the outer radius. If 6 is sufficiently small compared with a and if the condi- 
tions of support permit it it is generally considered allowable to obtain deflections and 
stresses in the ring plate by beam theory rather than by plate theory. 

The purpose of the present note is a more detailed analysis of this problem. We 
shall show that the question of whether beam theory and plate theory give essentially 
the same results depends not only on the ratio b/a but also on the magnitude of the 
deflections of the plate. 

Briefly, let ¢) be a quantity of the order of magnitude of the deflection of the plate 
according to beam theory divided by the width 6 of the plate, and let « be a dimension- 
less parameter of the form [12(1 — v*)]'”* b*/ah, where v is Poisson’s ratio and h the thick- 
ness of the plate. 

We find that beam theory is applicable as long as ¢u «K 1. When ¢. = O(1) then 
plate theory must be used and the appropriate equations are those of finite bending. 
We find further that when ¢ ou > 1 a boundary layer effect is encountered. 

In order to have a continuous transition between finite bending of beams and finite 
bending of plates we must use a system of equations for bending of plates which is more 
general than the equations of Kirchhoff and von Karman for small finite plate bending. 
\ more general system of this nature, applicable to axi-symmetrical bending of circular 
plates has recently been given.** 

The results which are obtained in what follows may be of some direct practical sig- 
nificance in connection with the analysis of expansion joints and of corrugated cy- 
lindrical shells. 

2. The equations of finite axi-symmetrical transverse bending of plates of constant 
thickness. The differential equations which must be solved are of the form** 


¢”’ + ! ¢’ — : cos ¢ sin ¢@ = ; [Vv sin jin (rV) cos ¢], (1) 
; ; " 


a ? ] ee 

yp’ 4 yy’ — ; cos @ — -¢ sn @¢ |V 
} r r 

(r' pu)’ 


— vpy COS d 
7 


0 
= — [cos¢@ — 1] — 
, 


a E cos @ sin @ + - ¢’ cos for) + ~ sin d(rV)’. (2) 
r r 


*Received Oct. 22, 1951. The present paper is a report on work done under the sponsorship of the 
Office of Naval Research under Contract N5-ori-07834 with Massachusetts Institute of Technology. 

**1). Reissner, On finite deflections of circular plates, Proc. Symp. Appl. Math. vol. 1 (1949) pp. 213-219 
and On axisymmetrical deformations of thin shells of revolution, Proc. Symp. Appl. Math. vol. 3 (1950) 


pp. 27-52. 
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In equations (1) and (2) r is the radial distance of a point of the middle surface of 
the plate before deformation, ¢ is the sloping angle of the deflected middle surface, ¥ 
is a stress function, p, the horizontal load intensity, V the vertical stress resultani, 
primes indicate differentiation with respect to r, and D and C are defined in the usual 
way by 
C = Eh, D = Eh*/12(1 — v*). (3) 
Stress resultants and couples and displacements of the points of the middle surface 


of the shell are given as follows 


rV = | rpy dr, rN, = VW cos¢ + (rV) sin 9, 
rQ = V sin 79) + (rV) COS Q, No = + pu 
‘ (4) 
M D(6’ “ * sin ), M,=- v{ sin @ + 16’) 
. r i 
u = - (NN, — vN,), w= i sin ¢ dr. 


The well known equations for small finite deflections follow from (1), (2) and (4) 
by writing sin @ ~ ¢, cos ¢ ~ 1 — }3¢, by omitting third and higher power terms in 
the dependent variables ¢ and WV and by omitting terms containing V in equation (2). 

3. Bending by transverse edge forces. In what follows we assume that 


Pu = Y, Py = 0. D) 

We write further 
rV = Pa, 6) 
so that P is the transverse edge load at the boundary r = a of the plate, per unit of 


circumferential length. 
We introduce dimensionless variables in (1) and (2) by writing 


r=ac+t bz, o = dof (2), VW = Wog(z). 7) 


This gives 


od : l 2,2 - Pa ( Baia ’ 
1 — of tee'Ifg- 1 — i coe], (8 
Dia + bea | 6 Po) fg Dia + bz) \ 2 PoJ + ) 
ei .. b ; b(1 — dof + ---) vbdi(l — sbof + ---) ff’ | 
—igqg’ + i — — — ( 
b } ' act ba 9 (a + bx) (a + bz) 9 


bof(1 — Fbof? + -+-) , vbof"(l — bof + --:) 
(a + ba)~ b(a + bx) 


7 _ 2.36 oe 
rs Vo! | l 6Po) - Pa, (9) 


(a + bx)’ 


where primes now indicate differentiation with respect to 2. 
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We now assume that the ratio b/a of width to inner radius of the plate is negligibly 
small compared to unity.* With this assumption equations (8) and (9) may be written 
in the form 


1 — Wob vr Ste 
/ ~ Da ( -* 6 Go] - ss ay eee oD \J pd 9 do) + °°], (10) 
| = | Ch o> ( A. ] 2 -2 9 ) “2 P boo ( = 1 2 2 : ) , 
Made l 12 dof + i + ." I 9 dof + “@ (11) 


Inspection of (10) and (11) indicates that the assumption b/a < 1 implies that the 
equations of beam theory can-be used for the solution of the ring-plate problem in the 
linear, small-deflection range. 

In order to see to what extent the linear small-deflection theory applies we may 
proceed as follows. We set 

Pb’ 
= ], (12) 
oD 
and we take account of the fact that.in non-linear plate theory non-linear terms both 
in (10) and (11) are important by setting 


V,b° ( b°do 


_ , 13 
Da V,a (13) 
From (12) and (13) follows 
Pb’ Eh? : 
— VY, = ——- 1 14 
“ D’ = lisa —»)"* (14) 
and 
Wb? 09/2 b? e 
—— = [1211 — »°)] - Ss li 
Da [121 v)I ah ? (15) 
Pbdy I h ‘ 
= 5172 - 16) 
V [12(1 —¥)]7 bo? 
We set as an abbreviation 
er sn 
= 9 — p _ 
iv ([12(1 v’)| a (17) 


and take account of the fact that always h/b < 1 and ¢ = O (1). Equations (10) and 
(11) then reduce to the following system 


fl’ - —(1 - Say" + ++) + ud — goof + fa, (18) 
9s I l jepe 2 9 
gh = gud — pods +o (19) 


*This is equivalent to saying that the leading terms only are retained in a development of the solution 


in powers of 6/a, 
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The following conclusions may be drawn from equations (18) and (19). 
(1) The condition fer beam theory to be applicable is of the form 


Udy K | (2) 
(la) When in addition 
oo <1 (72) 
then linear beam theory is applicable. 
(2) When 
ud, = O(1), (772) 


the governing equations are those of non-linear plate theory no matter how small 
ren) by itself may be 
(2a) Applicability of the Kirchhoff-von Karman theory of small finite deflections of 
plates requires satisfaction of (77) 
Let w, be a measure of the transverse deflection of the plate. We may set w) = @yb. 
Condition (i) becomes 
b l a , 
Wy < ao (120 — »)] y, Its v’) 
indicating that for the present plate problem linearity does not require that the de- 
flection be small compared with the thickness of the plate. 
4. Perturbation method for moderate values of ud, . In order to see for specific ex- 
amples to what extent the non-linearity for small ud, manifests itself the following series 


developments may be used 


20) 


| 


J = wbogi + (udbo) Js + 
In what follows we shall for simplicity’s sake assume that ¢5 by itself is negligibly small 
compared with unity so that the differential equations (18) and (19) may be taken in 
the form 


f"’ = —1 + udofg, g’’ = — judo . 21) 


We take as example a plate of width b with horizontal slope at z = +3 and with no 
horizontal edge forces. These boundary conditions correspond to what would be en- 
countered if the ring plate were part of a corrugated cylindrical shell. 

We have then as boundary conditions for the solutions (20) of (21) 


f.(+3) = 0 g(t 3) = 0 22) 
We find 
pelesg 
n= — h(t Set 4 Ba —H) - 
7 or a 7 ++ : .= a 2 a = a ~ im) 
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The relative edge deflection 6 is given by 


e1/2 


sa fs ¢¢= Be | fa 
- 20) | "fy dit + (wobe)* [- fe dx + +f (24) 
From (24) follows 
5 = te [1 — 0.000059(ud»)? + +--+]. (25) 


We see that while order of magnitude considerations require that ud) < 1 in order that 
there be no noticeable non-linear effect on load displacement relation the actual com- 
putations show that this is numerically conservative and that actually ud, can be ap- 
preciably larger than unity before an appreciable non-linear effect occurs. 

Having f, , fz and g, we may if we wish also determine the influence of small non- 
linearity on bending and direct stresses in the plate. We obtain in particular 


o, p(3) = 6M.(44)/0? = 4(3bP/h*)[1 + 0.00005(ud.)” + ---], 


oy pi +3) = No(+4)/hk = +(3bP/h*)(0.0012(1 — v*)**(udo) + ---]. 


5. Boundary layer equations for large values of ud, . In order to investigate the 
character of the solutions of equations (18) and (19) when yd, > 1 we introduce new 
dependent and independent variables by setting 


c=y, f=foF(y), 9 = g Gy). (26) 


We have then as differential equations 


F”’ = -[1 — : (bofo)’F? + | 3 + Xu - F (bof) PF? + ... |re, (27) 


2 ac. 
l ee > 
G 5A ud J [ —_ (dofo) + - | (28) 
2 g > Bi 
We set 
: oe . 
_* 1, N'udboJo = 1, N'udo  < l (29) 
oO! 
= (ud) *”* fo = (uo) *” Jo = (who) ”’. (30) 


Differential equations (27) and (28) now read 


Fr = {1 ~ iter 4 | + E wes ..- |r (31) 
Zu O pu 


I 


9 2 u 
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The order of magnitude of the change of slope ¢ is now ¢ofo = +/¢,/u and consequently 
we know that 

oo/u = O(1). (33) 
Accordingly all coefficients of (31) and (32) are O(1) and differentiation with respect to 
y does not change the order of magnitude of the quantities which are differentiated. 
Consequently significant changes of F and G occur over y-distances which are O(1). 
But this is the same as saying that significant changes of f and g occur over x-distances 
which are O(A) = O({udo]”'“*). Since udo > 1 we have then the existence of a boundary 
The actual width of the boundary layer with respect to the 7-coordinate is 


1/4 


layer.* 
b(udo) 

The question may be asked in which way the relative deflection 6 of the edges and 
bending and direct stresses depend on @» and u for the limiting case udp — ©. 

We may, again for simplicity’s sake, assume that ¢)/» is negligibly small compared 


to unity and restrict attention to the system 


F” = —1 + FG, G” = -—}F’. (34) 
If we take the same boundary value problem as in the preceding section but now let 
the edges of the plate be at x = 0 and x = 1, then we have at the inner edge 

e=y=0: F=0, G=0. (35) 


At x = 3 we have the symmetry conditions of vanishing F’ and G’ or with y as inde- 


pendent variable 
l 


a ert “6 OF’ = @ 
y= Dr = 2 Uo) : = U, 


G’ = 0. 36) 


As \ tends to zero the outer boundary y = 47‘ tends to infinity. We expect that for 


sufficiently large values of \~* we should be able to replace condition (36) by the limiting 
condition 


F’(o) = G’(o) = 0, (37) 


and the problem then is to solve the system (34) subject to the boundary conditions 


(35) and (37 
We expect that F’(y), G’(y) are O(1) and consequently bending and direct stresses 


have orders of magnitude determined as follows. 


ee h .,dofo Eh do », 
C-zR = 9 ko (r) = 9 y rb I \y) = 2b Pk I (y) 
or 

Oo, I ry" (2) F 

Dg - 9/ a= ght Ma 
E 2 (Z at "4 Nh ey 

(38) 
W'(r) — VoGo Wy) = « Eh 0 -Q"(y) 

ee ae ~ hb ry) = bf120. — a | asia a" ny 


*Boundary layer solutions for quite different problems of finite bending of circular plates have 
previously been discovered by Friedrichs and Stoker. See Am. J. Math. 63, 839-888 (1941). 
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ey = () f12(1 — v*)]' (2) G'(y) (39) 


It is remarkable that now o»p and o,, are of the same order of magnitude. Furthermore 
we see that in this range the stresses vary with the three-fourth power of the applied 
load and are independent of the actual width of the plate. They are inversely pro- 
portional to the thickness of the plate and proportional to the fourth root of the inner 
(or equally well average) radius of the plate. 
In order to determine the relative deflection of inner and outer edge we must calculate 
1/4 2A 


. »1/2A al 
6= | odr=2] — gofF(yrdbdy= 20%, | Fly) dy. (40) 
‘ /0 be 


Vow, before letting \ tend to zero we must know the behaviour of F(y) for large y. 
It seems that one might reason as follows. Assume that for large y 
F(y) = cy’, Giy) = cy’, (41) 
where both n and m are smaller than unity in order to satisfy (37). Equations (34) 
then become 
enn — 1)y = —1+¢eoy""”, cm(m — ly" * = —43ey. (42) 
The first of these equations gives 


C.Co = I, n+m=0, (43a) 


and the second c1ves 


c.m(m — 1) = —3c¢,, m — 2 = 2n. (43b) 


Equations 13) are satisfied by 


‘ 9 | 1/3 (*) 5 
m = = — C= », = ° }) 
3) n 3) Cc; (4) p ( (4 


We have then for large y 


A\* 273 9\'" 
F(y) ~ (4) ss G(y) ~ 1. (45) 


e 


If the first of these two formulas is introduced into equation (40) for 6 we conclude that 
since (45) should be valid over a predominant portion of the interval of integration we 


tbo 


should have 


; ol ry . | 1/3 , _¥ we. 4 ' 4] ( 1 ) 3 ™ (2) 3 ) 3 R 
6 ~ 2b " | (4) y dy = 2b e i 3(3 Dr = 6 9 b y 3 (46a) 
6 2\'"/ P\'“*/a\? a 

~ 6 461 
b a(?) (*) (2) — 


Clearly, the argument beginning with equations (41) is not rigorous. The results how- 


ever are not implausible, and one may hope that future more rigorous considerations 


will confirm them. 
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—NOTES— 


AN INEQUALITY FOR THE AMPLITUDES AND AREAS IN 
VIBRATION DIAGRAMS OF TIME-DEPENDENT FREQUENCY* 


By PHILIP HARTMAN ann AUREL WINTNER (The Johns Hopkins University) 


Let w(t) be any real-valued, continuous function on a ¢ interval, and let x(t) be any 
real-valued, non-trivial (40) solution of 


x’ +o(tzr = 0. (1) 


Unless w(t) is of some extremely particular type, no solution x(é) of this general form, 
|), of the problem of frequency modulation can be found “explicitly” (the actual case 
of frequency modulation is the particular case 


w(t) = 1/(a + B cos 0), (|B| <a); (1’) 


but (1) cannot be integrated in explicit terms in this case either). 

It is therefore of interest that, without solving (1), one can delimit an explicit in- 
equality for the maximum amplitude and the area of every half-wave (if any) of an 
arbitrary solution of every differential equation of the form (1). This is the content 
of (4) below. It will also be shown that the inequality in question is the best possible 
of its kind (if w(t) is unspecified). 

Unless the given /-interval is either too short or such that w(t) is “too close to 0” on 
it (i.e., (1) is “too close to x” = 0,” with x(t) = ¢,f + ce), a solution x(t) will have more 
than one zero, say x(a) = 0 and 2(b) = 0, where a < b. Then, since 2(é) is not the 
trivial solution, it can be assumed that ¢ = b is the first zero following the zero t = a; 
so that x(t) ¥ 0ifa <i <b. Since —2(t) isa solution of (1) if x(¢) is, there is no further 
loss of generality in assuming that 


a(t) > Oifa <t < b, while x(a) = 0, x(b) = 0. (2) 


Let M denote the maximum amplitude, and A the area, of such a half-wave in the 
(t, x)-plane, that is, put 


ab 
M \= max x(t), A= | x(t) dt. (3) 
a<t<b “a 
It will be shown that, no matter what the frequency function w(t) of (1) may be, the geo- 
metrical characteristics (3) of every half-wave (2) of every solution x(t) are subject to the 
inequality 


ab 


M/A <h | oo dt. (4) 


m= a 


It will also be shown that the 3 in (4) is the best absolute constant, in the sense that 
the } cannot be replaced by any numerical constant less than 3, valid for every differential 


equation (1). 


*Received October 10, 1951. 
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First, it is clear from (1) that 2’’(f) is non-positive if and only if x(é) is non-negative. 
It follows therefore from (2) that the (¢, x)-graph represents a convex arch over the 
segment [a, b] of the /-axis. Needless to say, the line ¢ = c, where c denotes the abscissa 
belonging to the maximum (for a < ¢ < 6) of the ordinate, need not be a line of :ym- 
metry of the arch. But it is clear from the convexity of the latter that the triangle 
having the base fa, 6] and the vertex ¢ = c, x = 2(c) is situated beneath the arch. In 
particular, this triangle has a smaller area than the convex region surrounded by the 
segment fa, b] and the arch. In view of (3), this means that 
A > 3(b — a)M. (5) 


On the other hand, by a well-known inequality which goes back to Liapounoff,' 
| w(t) dt > 4/(b — a). (6 


Clearly, (6) and (5) together imply (4) 

It is also clear that the inequality (5) turns into an inequality when the convex arch 
degenerates into two sides of the triangle considered above. On the other hand, the 4 
occurring in (6) is the best possible absolute constant,” and the optimal nature of this 4 
ean also be concluded by confining the convex arch to be arbitrarily close to the tri- 
angular form.’ Accordingly, the best absolute constants, } and 4, occurring in (5) and 
(6) respectively, can be approximated by the same sequence of examples. This proves 
the statement italicized after (4). 


iSee G. Borg, Amer. Jour. Math. 71, 67-70 (1949). For a proof embedding (6) into a general theory, 
ef. A. Wintner, ‘hid. 73, 368-380 (1951). For refinements of (6), ef. P. Hartman and A. Wintner, ibid. 73 
885-890 (1951 

E. R. van Kampen and A. Wintner, thid, 59, 270-274 (1937). 


’ 


A FORMULA FOR THE NORMALIZATION CONSTANT IN 
EIGEN VALUE PROBLEMS* 


$y GOTTFRIED GUDERLEY 
Given the differential equation 
y’ + Ar(a)y = 0 (1) 


with suitable homogeneous boundary conditions at the end of an interval which extends 
from « = ato x = b. In (1) \ denotes an eigen value. During the following derivation 
\ is considered as arbitrary and only the boundary condition at one end of the interval, 
say point a, will be imposed, so the functions y need not be eigen functions. If y = y; 


Received August 31, 1951. 
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and y = yr; are two such solutions for \ = A, and \ = Aj, respectively, then one has 
the relation 


2 Yuyr — YY : 
r(x)yy =* ———— 
| YY Ga he es es 


“a 


a 


(for a derivation see Ref. 1). Considering y as a function of z and 3, one obtains by 


the limiting process \;; — Ay , Yn — Yr , that 


3 r(a)y? dx = (24.3 = oy -y) , 

‘sf . Ox OX Ox Od * a 

If now the function y is an eigen function and its dependence upon d is explicitly known, 
this formula yields the normalization constant [; r(x)y’ dx immediately. Examples are 
the Jacobi polynomials or the functions of Ref. 2. In either case the functions y are 
expressed by hypergeometric series, the endpoints of the interval are singular points. 
The formula may even be useful in problems for which the eigen functions are deter- 
mined by numerical integration for different values of \ and subsequent interpolation. 


REFERENCES 
1. R. Courant and D. Hilbert, Methoden der mathematischen Physik, Interscience Publishers, New York, 
Vol. 1, Chap. V, See. 3. 


2. G. Guderley, Two-dimensional flow patterns with a free stream Mach number close to one, A. F. Technical 
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A NOTE ON THE HODOGRAPH TRANSFORMATION* 
By A. R. MANWELL (University College Swansea) 


The hodograph transformation [1, 2] for plane compressible flow is too well known 
to need any discussion here, whilst its singularities have been very fully investigated in 
(3, 4, 5]. Briefly, it has been found that in cases where plane potential flow is impossible, 
the continuation of the solution in the hodograph plane may be regular but cannot be 
mapped into the physical solution. In any given hodograph solution there is no difficulty 
in deciding if a solution is mappable, for the condition of the vanishing of the Jacobian 
of the transformation is just that the level lines of the stream function in hodograph 
space should touch the fixed characteristics in the hodograph space. This criterion of 
course requires the calculation of all the level lines. Beyond this not much is known, 
but it has also been shown by various writers, particularly Friedrichs [6], that limit 
lines cannot appear initially in a certain sense at points inside the solutions. The purpose 
of this note is to show that the somewhat abstract results of Friedrichs, may be included 
in the simpler and rather more general statement, that for a wide class of hodograph 
solutions, if the boundary streamline can be mapped then so can the whole field. 


*Received June 15, 1951. 
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1. The hodograph equations in characteristic form. The equations of isentropic plane 


potent ial flow are 


(pu) pv) 0 (1.3) 
a v 0 (1.2) 
dp l 
[ + 5(u +) = const., (1.3) 
J - Z 
dp 
(1.4) 
dp 
If 
u y p —=<_ (iS) 
U ¥./p Y 
and (1.6) 
u = qg cos @, v= q sin Gg. 
It is an elementary calculation to verify the equivalence of 
1 f (pq) \ ate 
op, + Pas y,> = 0, (1.7) 
J |" pq J 
l f q 1 
Lf v.? 0 (1.8) 
J \ he 
where 
1 , d= q C) 42 
J=p'aqivt ; Wo (1.9) 
(c.f. [4] See. 103). 
Another elementary calculation leads to the characteristic equations equivalent to 
(1.7), (1.8) viz. 
U. + 4B(V — VU) = 0, (1.10) 
V; + 3B(U — V) = 0, (1.11) 
where 
2da = dr/r — dd = q ‘(¢ c — 1)'” dq — dé, 
2d8 = dr/r + dé, 
and 
U = Wz h y = Wa /h 
where h*(q°/c’> — 1) = p’. 
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Also, 


and B is positive if it is agreed that the sign of the radical is to be chosen so that r in- 
creases with q. In this case both a and 8 increase with g. Now in the supersonic region 
the sign of J UV/q is that of UV, and for a continuous streamfunction defined in 
a simply connected region which is partly subsonic, the problem of showing that the 
solution can be mapped into physical space reduces to showing that UV is positive in 
the supersonic region. It will henceforth be further assumed that U and V are also 
continuous and piecewise analytic. Then if J changes sign it must vanish and, according 
to the theory of hyperbolic differential equations, discontinuities in the derivatives of 
(’ and V can occur only across characteristics, so that there will be Taylor expansions 
for U’, V in each quadrant of the a, 8 plane near the zeros of J. It will however be suffi- 
cient for the present purpose to consider the behaviour on the characteristics themselves. 

2. Expansions near a zero of J. Let a, 8 be, for the sake of convenience, measured 
from the point where J = 0. The cases that only one or both of U, V vanish must be 
considered separately. 

Suppose for example that U # 0 whilst V = 0 and along a = const = OV = 
V (8). Equation (1.11) shows that 


V(s) = —3B.U 8, (2.1) 
where the suffix 0 refers to the zero. Therefore 
and since B > 0 everywhere the sign of J, assumed positive in this paper for the sub- 
sonic region, is fixed with respect to increasing q. 
Similarly, if U 0 but V ~ 0, 
U(a) Via) & —4B Via. (2.3) 


If on the other hand both U and V vanish and so V(a) = V,a”" with n > 0, (1.10) gives 
U(a) > —4V,B,a"*’ so that 


U(a) Via) = —4ViBia"*’ (2.4) 
whilst similarly from (1.11) 
U(8) V(8) + —3V2B8°"*" 


It is perhaps necessary to point out that in the full Taylor expansions of U, V there 
may well be terms of lower order than m or n and that this is associated with cusped 
limit lines enclosing the axes a = 0 and 8 = O near the zero. Equations (2.2), (2.3), (2.4) 
show that J is strictly decreasing with increasing q along at least one of the characteristic 
directions through a point where J = 0. Its behaviour along general directions cannot 
however be predicted. 

At a sonic point the characteristic equations (1.10), (1.11) cannot be used directly 
but it is still desirable to retain this form as far as possible. To do this the limit of the 
previous equations will be considered for points near the sonic line. 

Let gq = c,(1 + é) where c, is the critical speed and ¢ small and positive. 
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Chen e ~e,{l —@ — Dé 
1/2 
h= p,/[(l + yt], 2.5) 
r=r(l + Ot )] 
B= HA + yey, 
while 
l 
U~py ive + (lh + ye] vo}, 
2.6 
F ar 1/2 
V~p iv: — (1 + weve}, 
and 
da = dg a l + y)t] dt. (2 7) 

If J vanishes (y 0 and near the zero on the sonic line the expansions along either 
characteristic are of the form U = Pt" and V & Qt" where m and n are not less than 

1.10 (or 1.11) after multiplication by @ and integration from 0 to ¢ gives either 
m= nandQ 1) 1)P or m > n and (4n — 1) 0; the second possibility must 
he rele al nN 

Pa 1? l prol Cc (2.8 
and there is a similar but in general distinct expansion on the other characteristic. 

Fre m 2 a 2 t 2.8 it appears that near a zero of e iL change 2 sign positive i” negalre 
along at least one characteristic in the direction of increasing q. The absolute sign in the 
subsonic region must of course be taken as positive. 

This contains the basie result found in [1] that near a zero of J it must change sign, 


and the conditions on the solution here are slightly less restrictive. 

Let the flow, or one sheet of it, be bounded by a smooth non-characteristic streamline 
C on which it will be further assumed that dY/dn is of one sign in the supersonic region. 
C will have one-one correspondence with a, 8 and J will be positive on C. Let a@,, be the 
largest lue of a for which J vanishes in the solution inside C and choose on a a 
the particular zero A with maximum 8. Now K cannot lie on C and from the form of 
the region in the a, 8 plane it is possible to join it to C by two characteristic ares each 
in the direction of increasing g. On these ares J > O by the construction whilst it has 
already been proved that near Bm a is essentially negative. It is therefore impossible 
for J to vanish inside the region. 

Figure l(a) shows the impossible case of isolated zeros disproved by Iriedrichs, 
Fig. 1(b) the more general situation excluded by this theorem 

There are two simple corollaries to this conclusion: 

i) For solutions of the above type the streamfunction cannot have a maximum or 
minimum at interior point s of the region. It is only necessary to show this for the super- 
sonic part and there if Y, = 0 and y, = 0 then J = 0 also which is impossible. 

ii) The uniquens ss theorem for Cauchy data on part of C. 

Let (U’, , V,;) and (U, , V2) be two solutions of the above type and if at some point 
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CHARACTERISTICS —~ 







BOUNDARY 
STREAMLINE 





l(a). Hodograph plane: J cannot be positive on the entire broken line in the vicinity of J = 0. 


CHARACTERISTICS 


BOUNDARY 
STREAMLINE - 


SONIC 
CIRCLE 





ig. 1(b). Hodograph plane: Limit lines, shown dotted, enclosing shaded regions of J < 0. 


inside the region Ul’, # U, (say) i.e. U, = AU, with A ¥ 1 the solution 


- = _[U,- dU, Vi - AV | 
w, 0) =| is Se eS 


has the same Cauchy data on the are C as the original solutions. However, by construc- 
tion, J for the new solution vanishes at one point inside the flow whilst the above 
theorem applies giving a contradiction. The region in which uniqueness has been proved 








182 NOTES [Vol. X, No. 2 


may be a triangle lying entirely in the supersonic region or it may be a quadrilateral 
extending from the are of C to the sonic line. (Fig. 2). 







BOUNDARY 
STREAMLINE 


Peeper eeel yf Pe Oe eee eee eee aes 


Fic. 2. Characteristic plane: Broken lines indicate typical boundaries for uniqueness with Cauchy 


data on the corresponding streamline arc. 


3. Solutions in the physical plane. If a hodograph solution has been found which 
satisfies the mapping condition then an infinity of neighbouring solutions can be con- 
structed which are also mappable. For, let ¥,(q, 6) be another solution having the same 
singularities and branch points in the subsonic region and, if ef, is added to the original 
solution, in general the boundary streamline and the values of the Jacobian on it will 
be altered, but, by quantities which by choosing e may be made as small as is desired. 
However, it was pointed out in [7] that infinitesimal changes of the physical boundary 
demand in general finite changes of the hodograph. This, which was shown in detail 
for subsonic flow is generally true. It is a consequence of there being no unit of length 
in the equations of motion of an inviscid fluid so that flow patterns do not depend on 
the scale of figure. Any small disturbance of the physical boundary is equivalent locally 
to the insertion of one half of a symmetrical aerofoil into the wall of a wind-tunnel. 
Since it is known that not all symmetrical aerofoils permit continuous flow it follows 
that general small displacements of the boundary will lead to shock waves. For sub- 
sonic flow however the variations must have large changes of slope and the disturbances 
will be merely local. 

In the case of local supersonic flow the situation is rather different. That some small 
variations of the physical boundary will lead to an end of continuous flow is a conse- 
quence of the existence of any symmetrical profile e.g. a very flat double wedge over 
which continuous flow has been proved impossible. Moreover the disturbances will be 
propagated along nearly characteristic paths in the physical plane to finite distances 
(Ref. 5, pp. 54, 131). Thus, whilst in the strict mathematical sense all compressible 
flows are ‘“‘unstable” in that near boundaries cannot have shock free flow, for transonic 
flow the “‘instability”’ has a physical meaning. C.f. [8, 9]. 

By means of the theorem of this note the situation for flows with simple hodographs 
xan be examined in more detail. (In physical flows there may also be both simple waves 
and flows needing several sheets in the hodograph). In any neighbouring solution corre- 
sponding to a small bulge in the physical plane only a few streamlines are affected. In 
the hodograph, which is supposed simple, these must initially pass into a new part of the 
plane. Also since the disturbance of the whole flow is small they leave through a small 
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gap in the hodograph boundary. On the other hand for finite changes of slope over the 
bulge they then must spread over a finite region of the hodograph plane. With these 
requirements it is, see Fig. 3, not possible to draw a smooth mappable hodograph 


boundary for the disturbed flow. 







CHARACTERISTICS 


BOUNDARY 
STREAMLINE 





SONIC 
CIRCLE 


Fic. 3. Hodograph plane: Limit lines appear at L,M on a typical perturbed boundary streamline, 
shown by broken line. 





Leaving aside the question of “stability” which, in situations in which shocks occur 
in a continuous manner, may not be of vital importance, one may enquire what theo- 
retical limitations the theorem of this paper puts on simple hodograph solutions. The 
answer to this question depends on the shape of the nose of the aerofoil. Usually this 
is normal to the stream and the boundary velocity up to the mid-section is turned 
through 90°. Comparison with the limiting characteristic shows that ¢/c* is limited to 
2.24. If the nose of the aerofoil is cusped inwards the boundary velocity can turn through 
180°, and the velocity increase found by tracing a characteristic through an angle 
(1/6 — 1)7/2 radians [130°] gives g/c* = +/6, the value for cavitation with y = 1.4. 
Therefore the theorem places no new limit on the local Mach numbers possible in tran- 
sonic flow and it seems quite likely that mathematical solutions with cavitation do in 
fact exist. 

Conclusion. The author first proved the theorem of this paper for the case that only 
one factor of J vanishes at each point on the limit line. This case, which is the usual 
one in applications, can be proved by showing that closed limit lines inside the flow 
are not possible whilst by hypothesis they cannot end on the boundary C. The proof 
given here was however arrived at after a study of reference [1] which also suggested 
the more general result given here. The discussion of the physical “stability” is largely 
independent of this theorem depending on the principle of the similitude of compressible 
flows at the same Mach numbers but different Reynold’s numbers viscosity being 


neglected in the mathematical theory. 
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NOTE ON A SUFFICIENT CONDITION FOR THE STABILITY OF GENERAL, 
PLANE PARALLEL FLOWS* 


By MARTIN LESSEN (The Pennsylvania State College) 


In reference 1, Synge derived sufficient conditions for the stability of plane Couette 
motion and plane Poiseuille motion. In the present note, it is demonstrated that similar 
conditions exist for any parallel flow with either finite or infinite boundary conditions 
or both. Although it can be shown that parallel two-dimensional motion must be of 
either Poiseuille or Couette type, it has been assumed, for purposes of stability con- 
siderations, that boundary-layer flows are substantially parallel in nature and Synge’s 
treatment has therefore been extended to cover these flows. 

The equation describing the stability of parallel flows for small disturbances is the 
Orr-Sommerfeld equation (reference 2): 

(w—c)(e’’ —ag) — wg = —- ; 7 (o'* — Qa’y’’ + a’y), (1) 
iL 


where z is the coordinate in direction of the flow, y the coordinate perpendicular to this 
the disturbance 


x(z-ct 


direction, ¢ the time, w = w(y) the velocity of steady flow, ¢(y)e’ 
stream function, c the disturbance phase velocity, a the disturbance wave number, R 
the Reynolds number of steady flow. 


*Received July 2, 1951. 
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If the boundaries (finite or infinite) of the flow field are at y, and y., the boundary 


conditions on the disturbance amplitude function ¢ are given by 


g(yi) = 0, g(y2) = 0, 


g’(y,) = 0, ¢g’ (yo) = 0. 


multiplied by ¢ dy, integrated over the flow field, 


If Eq. 1) is cleared of fractions, 
o its complex conjugate the following inequality 


and the resulting equation added t 
corresponding to Eq. (11.16) of Ref. 1 is obtained. 


alc, (1 at. a’ I;) = akg! yl, ep iz + 20° I? + a’ Is), (2) 
[- = | ge dy, i; = g’e’ dy, 
B= | oe" dy, Q =Winax c, = Ime. 


It should be noted that the values of all the integrals specified are convergent for 


infinite boundary conditions. This can be seen by considering the asymptotic properties 


of gat yo a in the integrals. 
o(y ~+«~) = constant X (e°“”) 
For a discussion of the asymptotic properties of ¢, see reference (3). 
; ¢ 


ving inequality is now considered in place of (11.17) of Ref. 1 because the 


latter cannot be used for boundary conditions at infinity. 
| |et+éie’ +n” |? dy >0 


ive real, arbitrary constants. It follows that 
nil; >(2n-®), -—Ih. (3) 


3) is substituted inte the inequality (2), the following relation results. 


Ren (I? — oe I) < aRqr Tol, — Ti(2n — & + 20°) — Ii(a'n — 1). 
(alqn’)” > 4(2n —E + 2a’ n)(a'n’ — 1), (4) 


where 


2n —& + 2an > O, a'n —1>0. 


The sufficient condition derived for stability of general, plane parallel flows is a very 
weak one but it has the advantage of being extremely simple. The condition derived is 
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weaker than the one obtained by Synge but it has wider application. Such two-dimen- 
sional, parallel flows as boundary layers, jets, and wakes can be quickly considered 
preliminary to more detailed investigation. 

Consider Eq. (4) as applied to the laminar boundary layer between parallel streams 
(Ref. 3). 


Let » = K/a’, # = 0. Equation (6) then becomes 


If the right hand side is maximized with respect to K, the corresponding value of K is 
K = 3. The value of g for the problem under consideration is 


q = 0.2000. 


Therefore, for certain stability 
R < 15.8a°. 


It is interesting to note that, although no lower branch of the stability curve is 
obtained by the above method, the actual calculations to a second approximation for 
the free boundary layer also do not vield a lower branch of the neutral curve. 
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GENERALIZATION OF A PROBLEM OF RAYLEIGH* 


By STANLEY CORRSIN (Department of Aeronautics, The Johns Hopkins University 


Emmons has recently published a note [1] giving dissipative temperature distribu- 
tions for Rayleigh’s classic problem of the incompressible viscous flow set up by an 
infinite plate suddenly moved at constant velocity in its own plane [2]. This suggests 
sufficient interest to warrant exposition of a particular solution of the problem for a 
compressible fluid. This was worked out several years ago as a simple pedagogical 
demonstration of some characteristics of the compressible boundary layer, much as the 
Rayleigh case is often used to introduce students to the incompressible boundary layer. 

We shall restrict the problem by taking low Mach number, v/u < 1, insulated wall, 
Prandtl number unity, perfect gas, u ~ \ ~ T°”* and c, = constant. The plate lies on 
the plane y = 0 and moves in the z-direction; u is velocity along 2, v is along y; ¢, , u, 
, p, 7’ are specific heat at constant pressure, viscosity, thermal conductivity, density and 
absolute temperature, respectively. 


*Received July 5, 1951. 
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The z momentum equation with v/u < 1 is 
Ou te) au) 
oad | Sed & l 
P at oy (. Oy (1) 
For low Mach number, the energy equation is 
aT a ( af) (2) 
eo hh ae 2 
> “at Oy Oy +* oy (2) 
and the state equation is 
p T; 
= 7 (3) 
Pi 1 
where the subscript 1 refers to y = ~. Also, 


(4) 
Assumption of unity Prandtl number permits satisfying (1) and (2) simultaneously 
with ¢,T = f(u) only [3]. Substitution of this into (2) shows that (2) is identical with 
(1) only if 


af 


du? 


= -1 


whence, 


¢,T = A+Bu-—<- (5) 
2 
The boundary conditions for u and T are u 0 and T = T, at y =@ ort = 0; 
i= UandT = T, at y = Oort =o. T, is stagnation temperature for 
with velocity U’ and temperature 7’, 
Then (5 


a fluid flowing 
since the plate is insulated. 
can be transformed to 


ee - i u \? 
T= T —(T) —7)(1—%). 





(6) 
Since there is no heat transfer between wall and fluid at y = ©, the one-dimensional 
energy equation is applicable: 
(iP 7 — 1 2 -_ 
T, L+“— 54h, (7) 
where the effect of pressure (~ dp/dt) has been neglected with the restriction to low 
Mach number, @, = U/a, and a, is the velocity of sound at y = ~; y is the ratio of 
specific heats. 
With (7) and the new velocity variable w = (1 — u/U), (6) becomes, 
i 


— 1 i 2 
771+ a t Mi(1 — @’) 


(8) 
which gives 7'(y, ¢) after u(y, t) is known. 
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With (4) and (8), the x momentum equation becomes 


Ww 5 ara CE 3 - ‘. ; Ow \” 
vy fl + NO — w)]* = — =v,N[1 + NU — &)] AY, (9) 
2 OY 


al oy 


where 


As in the incompressible case, this is reducible to a total differential equation with a 


similarity variable, 


j 
i 2(y,t 
Then 
daw 3 Nw du) 27 da 
} a - : i T 2 = } i: ! = = () 10) 
an 2 1+ N(1 —w) \dn |] a N(1 w )|' , dn 
with boundary conditions w = 0 at 7 O;w9 = latyn =m. 
Consistent with the restriction to small Mach number, we assume N(1 — “So 


and (10) reduces to a relatively simple non-linear total differential equation: 


da’ . , dw 0 | | 
als = ; ) 
dn 


This is easily integrated twice to give, after application of boundary conditions, 





erf [(3N)' “ew eri [((2N) “| ert (7). 12) 
For , — 0, this reduces to the well known solution of Rayleigh, 
= f 13 
- = el nH). o) 
{ 7 


The skin friction coefficient c, can be obtained directly from (10), without going to 


the simpler form, (11). The wall shearing stress 1s 
4) U (usp,\"" ; les 
. arr fj = (mes) (1 + N) (%) (14) 
OY o 2 fs dy/ 


To get (dw/dn), , we integrate (10) once and let » — 0. This gives 


or 
( at, = (5) (1 + 7 = Mi) 16) 


which shows a Mach number behaviour much like that of the boundary layer. 
It seems physically clear that when the plate starts, a pressure wave is generated at 
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the wall and moves outward. Professor L. Howarth has mentioned to me that he is 
working on this more complex aspect of the problem. [Ref. added Oct. 1, 1951: L. 
Howarth, Some aspects of Rayleigh’s problem for a compressible fluid, Q. Jour. Mech. 
and Appl. Math., 4, pt. 2, 157 (1951)]. 
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THE QUARTER-INFINITE WING OSCILLATING AT SUPERSONIC SPEEDS* 
By THEODORE R. GOODMAN (Cornell Aeronautical Laboratory, Inc.) 


Introduction. In the April, 1951 issue of the Quarterly of Applied Mathematics, two 
papers appeared. One, by H. J. Stewart and Ting-Yi Li’ purports to extend the method 
of J. C. Evvard’ for solving supersonic wing tips to the case of the oscillating wing. 
The other by J. W. Miles” solves the problem of the oscillating rectangular airfoil by 
the Wiener-Hopf technique. In an application of their method to the rectangular airfoil,‘ 
Stewart and Li obtain results which disagree with those of Miles. Moreover, Stewartson® 
has treated the problem independently, and his results agree with Miles’. 

The raison d’étre of this note is to add yet another independent solution to the list 
of those which agree with Stewartson’s and Miles’. In view of the fluxive state which 
this problem is in at the moment, it is felt that this would be of current interest. 


The Gardner method. The method presented by C. 8. Gardner’ will be used to obtain 


the solution. Let x, y, z represent rectilinear space coordinates, x in the direction of the 
free stream, and y in the “spanwise”’ direction. Denote by U’ the freestream velocity, c 
the velocity of sound, JJ the Mach number and / the time. Let z = f(x, y, 0) describe 
the surface of the wing, g = US, + f,,8 = (M? — 1)'?, 2’ = 2/8, = (Mx — p'et)/8. 
Then in terms of 2’, y, 2, “, Gardner’s method for the rectangular airfoil consists in 

Received July 25, 1951. 
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solving the following two boundary-value problems in succession: (The primes on 2’ 
/ ~ \ 
and ?’ have been disregarded). 


I) For &¢>0:x,, Xer — Xee = O. (1) 
For & = 0: x a(x, y, t). (2 
For zx < 0:x 0. (3) 
II) ¥ . . 0 (4) 
For y > 0: lim vy, = lim, VY, = x(E, 2, ¥, 2). (5) 
For &>a:y 0. (6) 

The velocity potential, y, is determined by setting & = 0 in the solution y. 


The first problem. The value of g is taken to be independent of y, and in terms of 
physical « — ¢ variables is given by 


g = U(x) exp (tw). (7) 
In terms of the transformed variables, this becomes 
gq = vo(x) exp [w(Mz — 0)], (8) 


where vy = wM/BU. 











Fic. 1. “Wing” Problem in x — ¢ Space. 


Let x = o;: ; then o may be determined by integrating sources over the shaded region 
shown in Fig. 1. 


DP ae 
o(t, x,t) = —— | — 09(2,) exp (Mz,) dz, 


exp (—ivt)) dt 





<7. (9) 


Je-teez,)2-e2p270 [(2 — 4) — (t — 4) -— & 
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If the substitution ¢ — t; = [(~ — 2,)” — €]'” cos 6 is made, the inner integral is 
recognized as a Bessel function. Then after some manipulation, Eq. (9) becomes: 


o(€, x, t) = —exp [i(Maz — 0] [ v(x — X) exp (—ivMX)Jo[o(X* — &)'?] dX. (10) 


vf 
‘ 


. . 7 ° ° 
The second problem. The author has obtained elsewhere’ some general relationships 
which arise in using Gardner’s method for the rectangular airfoil. In the two-dimensional 


region 
g = —a(0, 2, 2). (11) 
Hence, 
yg = exp [(Mz — 0] [ v(x — X) exp (—vMX)J,(vX) dX. (12) 
In the tip region 
eg = —00,2,) +2 'y'” | a(t, x, thé — y)'”” dé/é. (13) 


Hence, 
re exp [tv(Ma — 2)] | v(x — X) exp (—wMX)J(vX) dX 


— ry” exp [in(Maz — 0] | (¢ — y)'”* dé/é 


| vole — X) exp (—i»MX)Jop(X* — #7] dX. (14) 


The wing and regions are depicted in Fig. 2. 





TWO DIMENSIONAL REGION 


TIP 
REGION 
So MACH LINE 





Fic. 2. Diagram of Wing 


7Goodman, Theodore R., Aerodynamics of a supersonic rectangular wing striking a sharp-edged gust, 
J. Aero. Sci. 18, 519-526 (1951). 
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NOTES 


The spanwise integral of the velocity potential. Let y represent the integral of ¢ in the 


direction from the tip to the semispan s: 


Y) eCXp wllNX)JS,(vX) dX 


dy | (£ 


t)| [ y 


vMNX)J,[v(X° gai Fs 


mly the triple integral appearing in the last term. Denote this by I. Upon 
w the order of integration, there is obtained 
| Y) exp (—7WATNX)J [v(x | dX / ly /(é y) | dy 16 
( ing with respeet to y, there obtained 
t=x/2 | de | ole — X) exp (—ivMX)Jolo(X? — &)""] aX 7 
he ler of integration once again 
/ 2 | ve — X) exp (—ivMX) dX [ Fob(X? — #)*] db. (18 
x 6, the last integral will be recognized as a Struve funetion 
/ Da Y) exp ivMNX)H _,,.(vX)(r/2vNX)'” X dX 19 
A() 


H] vi 2 /mv.X sin v.\ 


hie ( btained Lol 4 
| Y) exp vi X)| J ,(vX sin v.X)/2sy] dX Z 
en for the difference in notation, it will be seen that this solution 
\l les’ nd Stew irtson’s 
| ly \L[nem Ilan Company, 1945, p. 374, I 
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BOOK REVIEWS 


Partial differential equations in physics. By Arnold Sommerfeld. Academic Press, Inc., 
New York, 1949 


I} book is essentially concerned with the treatment of the classical linear equations of mathe- 
matical physics. It opens with the expansion of a function in a trigonometric series and the extension of 
h tion to the Fourier Integral and generalized Fourier Series. The second chapter distinguishes the 
elliptic, hyperbolic, and parabolic partial differential equations and discusses the analytic character of 
the solutions of such equations. It also introduces the Green’s function and Riemann function and the 
ted methods of integration. Chapter III is confined to heat conduction problems and Chapter IV 

treats a large variety of boundary value problems having spherical or cylindrical symmetry. 
Chapter V deals with various eigenvalue problems and the final chapter is concerned with electro- 
gnetie radiation. The book contains an excellent choice of very instructive problems and should be 

invaluable in teaching the technique of solving the classical type boundary value problems. 


G. F. CARRIER 


Differential equations. By H. B. Phillips. John Wiley & Sons, Inc., New York, and 
Chapman «& Hall, Ltd., London, 1951. viii + 149 pp. $3.00. 


This is a text-book on the theory and application of the usual types of ordinary differential equations 


lementary solutions. Types included are: first order with variables separable, other first order 

type pecial types of second order equations, and equations with constant coefficients. Both methods 

ution and development of the equations for physical systems are treated. Applications are taken 

eral fields with emphasis on mechanics. The discussion of variable mass problems is particularly 

r. This book seems particularly suited for a first course in differential equations for science and 
engmeering student 


E. H. Lee 


Lectures on classical differential geometry. By Dirk J. Struik. Addison-Wesley Press, Inc., 


Cambridge 42, Mass., 1950. viii + 221 pp. $6.00. 


Chis text, written by an outstanding geometer, contains an excellent vector treatment of the classical 
theory of Metric Differential Geometry in Euclidean three-space. The author’s style is clear and concise 
ind offers the reader a deep geometric insight into the subject. In addition, the text contains: (1) a very 
good collection of problems and numerous illustrative examples; (2) more than one hundred excellent 
figures; (3) many interesting historical notes. 

In the first section of the text, the author discusses the theory of space curves and considers such 
topies as: curvature and torsion; contact; the Frenet-Serret formulas; the fundamental theorem; some 
results on ovals in differential geometry in the large. Secondly, the imbedding theory of surfaces in three 
pace is discussed. Here, the stress is on curvature properties such as: Euler’s theorem; the Dupin indica- 
rix, asymptotic directions; the fundamental theorem of surface theory. Next, the intrinsic properties of 
urfaces (geodesic curvature, geodesic coordinates, etc.) are considered. Finally, the following special 
topics are briefly discussed: conformal maps, minimal and ruled surfaces. 

N. CoBurn 
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Integral transforms in mathematical physics. By C. J. Tranter. Methuen & Co., Ltd., 
London, and John Wiley & Sons, Inc., New York, 1951. ix + 118 pp. $1.50. 


This member of the Methuen Monograph Series contains a concise account of the application to 
physical problems of the Laplace, Fourier, Hankel, and Mellin transforms. After defining the transforms 
and their inversion formuli, the use of each is illustrated by relatively conventional example problems. 
The techniques of finding asymptotic and power series developments of solutions defined by not easily 
inverted transforms are discussed and Filon’s numerical inversion technique is presented. 

Finally, the author presents a chapter on the combined use of the finite transform and the relaxation 


technique 
G. F. CARRIER 


An introduction to electron optics. By L. Jacob. Methuen & Co., Ltd., London, and John 
Wiley & Sons, Inc., New York, 1951. x + 150 pp. $2.00. 


The author has compressed a great deal of material, including a large number of diagrams, into this 
brief survey of what is now a subject of considerable magnitude, both on the experimental and theoretical 
sides. The book contains thirteen chapters, list of general reference books, bibliography chapter by 
chapter, and index. The subject matter falls into two roughly equal parts, the first eight chapters dealing 
with optical concepts and analogies, electrostatic lenses (their image-forming properties and aberrations), 
and phase-focussing, while the last five chapters deal with magnetic lenses and the deflection and proper- 
ties of beams 

The purpose of the series, as set out on the jacket, is “to supply readers of average scientific attain- 
ment with a compact statement of the modern position in each subject.” Such a hypothetical reader 
should undoubtedly get from the book a good general impression of what the problems of electron optics 
are, and much information of a more detailed nature, which he can supplement by looking up the papers 
cited in the bibliography. Nevertheless, the proper achievement of the worthy purpose of the series 
demands (in the opinion of the reviewer) much more forethought in arrangement and care in detail than 
more in fact than are necessarily demanded in a book intended for special- 
rt 


are to be found in this book 
ists, since the specialist has greater powers of digestion and correction than the hypothetical reader “ 


t 


average scientific attainment.’ 
It is safe to say that such a reader gets nothing but confusion from formulae in which the notation 

is insufficiently explained or which are marred by a general carelessness, and his plight is still worse when 

this carelessness is also present in the descriptive writing. Here are some examples: 

p. 4: “where s; , s2 is the length of path” should read ‘‘where s; , s2 , etc. are the lengths of the paths. 


p. 14: Eds is used for a scalar product, without explanation. 

p. 15: A minus sign is omitted at least five times in equations (IV.4). 

pp. 15, 16: The unsophisticated reader will be puzzled to read that the cur! of E is an important concept 
associated with the field, and a few lines further that E possesses no curl. 

pp. 16, 17: The usual symbol for partial differentiation is used unnecessarily for ordinary derivatives, 
and changed to d suddenly. 

pp. 37-39: The 3-page chapter on the trajectory equation deals only with motion in a plane of symmetry; 
the reader will be in a state of doubt as to whether this is the general] case. 

p. 58: The sentence ‘‘the beam angle has begun to exceed the value for which the sine condition sin 6 = 6 
holds” would be vastly improved by the addition of the word ‘“‘approximately’”’; but even then the 
words ‘‘the sine condition” make one wonder what the author had in mind. 

These things are in a sense very trivial, but not so for a reader of average scientific attainment who 
is struggling to understand a difficult subject. Like other subjects in mathematical physics, electron 
optics has a wide range, from description of experiments on the one hand to mathematical theory on the 
other, and no author is to be blamed for failing to do justice to the whole range in a single short book. 
But the very difficulty (or impossibility) of the task is all the more reason for avoiding gaucheries of the 
tvpe indicated above 

J. L. SYNGE 
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The classical theory of fields. By L. Landau and E. Lifshitz. Translated from the Russian 
by Morton Hamermesh. Addison-Wesley Press, Inc., Cambridge 42, Mass., 1951. 
ix + 354 pp. $7.50. 


This book constitutes an attempt to treat the theory of electromagnetic and gravitational fields in a 
systematic manner. The authors restrict themselves to classical fields; no mention is made of the quantum 
theory of fields, although there is a reference to the fact that a limitation is placed on certain aspects of 
classical electrodynamics by quantum effects. 

By way of an introduction to the theory of the electromagnetic field (which is presented principally 
in the four-dimensional notation of the special theory of relativity), the first two chapters are devoted to 
the principle of relativity and relativistic mechanics. These two chapters are written very clearly, but 
present only a review of those parts of special relativity which are required for the remainder of the book; 
they are far too condensed for a student not already familiar with the physical concepts underlying special 
relativity theory. 

The following seven chapters deal with classical electromagnetic field theory. The topics considered 
in the first two of these chapters include the motion of charged particles in uniform fields, transformations 
and invariants of the fields, Maxwell’s equations and the energy-momentum tensor of the field. Chapter 5 
includes a discussion of constant fields and dipole, multipole and magnetic moments. It also touches on 
the field of a moving charge, but this is rather the subject matter of Chapter 8. Chapter 7, entitled 
propagation of light, deals with various aspects of geometrical optics and diffraction. Finally, as far as 
electromagnetic field theory is concerned, Chapters 6 and 9 deal with electromagnetic waves and radia- 
tion of such waves respectively. 

The theory of electromagnetic fields is limited to vacuum and point charges, the problems of continu- 
ous media being completely omitted. It is given in both three and four-dimensional notation, with 
particular emphasis on the latter when there are definite advantages to be gained from its use. 

In the eyes of the reviewer there is one regrettable feature in this part of the book, namely, the purely 
formal manner in which the theory is developed with its consequent neglect of the physics. 

The last two chapters of the book treat particles in a gravitational field and the gravitational field 
equations. While an attempt is made to discuss the concepts underlying general relativity theory before 
launching into the mathematical formalism, it is the opinion of the reviewer that the authors pass off the 
subtle questions of the meaning of space and time in general relativity much too glibly. 

There is a definite need f-r a book that covers electromagnetic and gravitational fields and points 
out the similarities and differences of the problems encountered in each case. This book is a good step in 
this direction for the reader who is at home with the basic concepts and physical principles of electro- 
dynamics and special and general relativity. 

L. C. Maxton 


The algebra of vectors and matrices. By Thomas L. Wade. Addison-Wesley Press, Inc., 
Cambridge, 1951. ix + 189 pp. $4.50. 


This is a clearly written text, on an elementary-intermediate level, dealing with the algebra of 
vectors and matrices. In particular, the author’s treatment of matrix algebra is illustrated by many 
examples and should provide a good introduction to this subject. Although the point of view of the 
author is that of the mathematician, there are a sufficient number of interesting results in matrix theory 
to make this a useful book for the engineer. 

The contents of the book are divided into two sections: Vector Algebra and Matrix Algebra. In the 
first section, the author treats the scalar and vector product and linear dependence of vectors in two, 
three and n-dimensional Euclidean space. The section on Matrix Algebra covers: addition and multipli- 
cation of matrices; groups of linear transformations; the characteristic equation and the Cayley-Hamilton 
theorem; the rank of a matrix and linear dependence; reductions of matrices, and of bilinear and quad- 


ratic forms to canonical form. 
N. Copurn 
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Universal mechanics and Hamilton’s quaternions. By Otto F. Fischer. The Axion Institute, 
Stockholm, 1951. vi + 356 pp. $10.00. 


The quaternion retired from active participation in mathematical physics over fifty years ago, 
withered by the caustic wit of Oliver Heaviside (Electromagnetic Theory, Chap. III), and its reappear- 
ances since have been few and far between. As the title of this book implies, the author believes in the 
power of the quaternion. His Preface begins as follows: 

“This is a book written by a civil engineer on universal mechanics with an attempt to introduce a 
certain order in its mathematical structure by means of the calculus of Hamilton’s Quaternions.” 

“The term ‘universal mechanics’ refers to the mathematics of ordinary physics of motions, elasticity, 
hydrodynamics, aerodynamics, electromagnetism, together with relativistic and cosmic physics as well 
as quantum mechanics. With Eddington we may divide it into three main groups of atomic physics, 
ordinary ‘molar’ physics and cosmic physics.” 

“In attacking a subject of such magnitude with a new tool it is difficult to follow the regular roads, 
I have called the book a ‘cavalcade’ indicating a rather erratic path of approach. To begin with, the 
quaternions are applied in a rather haphazard manner to a series of examples, a wandering together with 
the reader through many ‘countries’ in which I have found new beauties and exhilarating outlooks. 
Gradually, however, new lines of structure begin to take shape, growing more and more clear until con- 
tours of mighty trunk lines seem to appear.” 

The book contains four chapters with the following headings: I. Extracts of ordinary vector mechan- 
ics. II. Simple quaternions in mechanics. III. Quadric quaternions and Eddington’s E-numbers. IV. 
Further developments and ideas. Chapter II takes up more than half the book, and some idea of its 
variety will be given by the following samples from its 66 headings: Quaternation or rigid rotation. Max- 
well’s equations and Minkowski’s space. Electromechanical parallelisms. Mobius’ null-system and 
Pliicker’s line coordinates. On the theory of the general magneton. Eddington’s electromagnetic energy- 
momentum tensor. On gravitation and restricted relativity. Einstein’s theory of gravitation. Einstein’s 
theory of 1950. 

The book is rotaprinted. This does not seriously detract from the legibility of the letterpress, although 
it is a little too small for full reading comfort. But the formulae are cramped, being written in by hand on 
the original typescript in spaces inadequate to receive them. The practice of marking vectors with a 
straight underline and quaternions with a wavy underline is good, but the placing of a suffix directly 
under an underlined symbol is objectionable because it gives the appearance of a fraction. One particu- 
larly unfortunate slip must be mentioned; in the formula (2.1) on p. 39, defining a quaternion, the 
quaternion is equated to a vector. 

This book will not, I fear, fulfil the purpose of bringing the quaternion into common use. For the 
merit of the quaternion lies in its algebraic simplicity, and a direct clear statement of fundamental 
properties, such as that given in L. Brand’s Vector and Tensor Analysis (Wiley, New York, 1947), is 
better propaganda for the quaternion than pages of complicated formulae in which the basic simplicities 
are lost. However, one should not judge a book in terms of its propaganda value—that is a reaction 
induced in the reader by over-enthusiasm on the part of the author and by a hangover from the old days 
of quaternionic controversy. The book is not, as the author remarks, one to be read straight through. 
It contains a mine of formulae, ranging over a wide extent of physics, and it is in fact a book to keep on 
one’s shelves and consult in connection with special interests. There is an adequate index and a list of 


references. 
J. L. SyNcE 


Table of the reciprocal of the gamma function for complex argument. By J. P. Stanley and 
M. V. Wilkes. Computation Centre, University of Toronto, 1950. 101 pp. $4.50. 


The reciprocal of the Gamma function I'(Z) is calculated in the rectangular domain —.5 < ReZ <. 
0 < ImZ <1, in increments of .01. The functional values are given to seven decimal places. 


G. F. CARRIER 
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